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A coarse-grained simulation model eliminates microscopic degrees of freedom and represents 
a polymer by a simplified structure. A priori, two classes of coarse-grained models may be 
distinguished: those which are designed for a specific polymer and reflect the underlying atom- 
istic details to some extent, and those which retain only the inost basic features of a polymer 
chain (chain connectivity, short-range excluded-volume interactions, etc.). In this article we 
mainly focus on the second class of generic polyiner models, while the first class of specific 
coarse-grained models is only touched upon briefly. Generic models are suited to explore gen- 
eral and universal properties of polymer systems, which occur particularly in the limit of long 
chains. The simulation of long chains represents a challenging problem due to the large relax- 
ation times involved. We present some of the Monte Carlo approaches contrived to cope with 
this problem. More specifically, our review contains two inain parts. One part (Sec. 5) deals 
with local and non-local updates of a polymer While local moves allow to extract information 
on the physical polymer dynainics from Monte Carlo simulations, the chief aiin of non-local 
moves is to accelerate the relaxation of the polymers. We discuss soine exainples for such 
non-local moves: the slithering-snake algorithm, the pivot algorithm, and its recently suggested 
variant, the double-pivot algorithm, which is particularly suited for the simulations of concen- 
trated polymer solutions or melts. The second part (Sec. 6) focuses on modern Monte Carlo 
methods that were inspired by the Rosenbluth-Rosenbluth algorithm proposed in the 1950s to 
simulate self-avoiding walks. The modern variants discussed comprise the configuration-bias 
Monte Carlo method, its recent extension, the recoil-growth algorithm, and the pruned-enriched 
Rosenbluth inethod, an algorithm particularly adapted to the simulation of atfi'actively interact- 
ing polymers. 



1 Introduction 

Polymers are macromolecules in which N monomeric repeat units are connected to form 
long chains." Experimentally the chain length N is large, typically 10"^ ^ N < 10^. 
The size of a chain lO^A) thus exceeds that of a monomer (~ lA) by several orders of 
magnitude. However, contrary to granular materials,^ the chain is not so large that thermal 
energy would be unimportant. Not at all! Thermal energy is the important energy scale 
for polymers. It provokes conformational transitions so that the polymer can assume a 
multitude of different configurations at ambient conditions.'^ 



"More precisely, this definition refers to "linear homopolymers", i.e., linear chain molecules consisting of one 
monoiner species only. By contrast, polymer chemistry can nowadays synthesize various other topologies, such 
as cyclic, star- or H-polymers. For a very commendable review on the physical chemistry of polymers see Ref. 1. 
''Here, we mean the thermal energy supplied at ambient temperature, i.e., k^T = 4.1 ■ 10~^^ J for T = 300 K. 
'^Polymers are a paradigm for "soft matter" materials or "complex fluids". Roughly speaking, "soft matter" 
consists of materials whose constituents have a mesoscopic size (microscopic scale ~ lA <C mesoscopic object 
~ 10^ - 10*A < macroscopic scale ~ 1mm) and for which k^T is the important energy scale (whence the 
softness at ambient conditions). Examples other than polymers are colloidal suspensions, liquid crystals, or fluid 
membranes.^ 
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Changes of the configurations occur on very different scales, ranging from the local 
scale of a bond to the global scale of the chain.^-^ This separation of length scales entails 
simplifications and difficulties. Simplifications arise on large scales where the chain ex- 
hibits universal behavior. That is, properties which are independent of chemical details.*'^ 
These properties may be studied by simplified, "coarse-grained" models, e.g. via com- 
puter simulations. For simulations the large-scale properties, however, also give rise to a 
principal difficulty. Long relaxation times are associated with large chain lengths.^"' 

The present chapter focuses on some of the Monte Carlo approaches to cope with this 
difficulty. Why Monte Carlo? Within a computational framework it appears natural to 
address dynamical problems via the techniques of Molecular Dynamics (see Ref. 10). A 
Molecular Dynamics (MD) simulation numerically integrates the equations of motion of 
the (polymer) system, and thereby replicates, authentically, its (classical) dynamics. As the 
polymer dynamics ranges from the (fast) local motion of the monomers to (slow) large- 
scale rearrangements of a chain, there is a large spread in time scales. The authenticity 
of MD thus carries a price: Efficient equilibration and sampling of equilibrium properties 
becomes very tedious -sometimes even impossible- for long chains. At that point, Monte 
Carlo simulations may provide an alternative. Monte Carlo moves are not bound to be 
local. They can be tailored to alter large portions of a chain, thereby promising efficient 
equilibration. The discussion of such moves is one of the gists of this review. 

Outline and Prerequisites. The plan of the chapter is as follows: We begin by gathering 
necessary background information, both as to polymer physics (Sec. 2) and as to the Monte 
Carlo method (Sec. 3). Then, we present the simulation models (Sec. 4), which have been 
used to develop and to study various Monte Carlo algorithms. The discussion of these algo- 
rithms (Sees. 5 and 6) represents the core of the chapter. Section 5 deals with local moves, 
allowing to study the physical polymer dynamics via Monte Carlo, and non-local moves 
(slithering-snake algorithm, pivot algorithm, double-pivot algorithm), aiming at speeding 
up the relaxation of the chains. Section 5 discusses the Rosenbluth-Rosenbluth method for 
simulating self-avoiding walks and some of its modem variants (pruned-enriched Rosen- 
bluth method, configuration-bias Monte Carlo, recoil-growth algorithm). The last section 
(Sec. 7) briefly recapitulates the different methods and gives some advice when to em- 
ploy which algorithm. Finally, the appendix A reviews a recently proposed approach to 
systematically derive coarse-grained models for specific polymers. 
Our presentation is based on the following prerequisites: 

• We will restrict our attention to homopolymers, i.e., to polymers consisting of one 
monomer species only. However, (some of) the algorithms discussed may also be 
appUed e.g. to polymer blends or block-copolymers (see Ref. 11). 

• The chains are monodisperse, i.e, N is constant. 

• We do not consider long-range (e.g., electrostatic) or specific (e.g., H-bonds) inter- 
actions between the monomers. These interactions are treated in other chapters (e.g., 
see Refs. 12, 13). 

• We do not treat the solvent molecules explicitly. They are indirectly accounted for 
by the interactions between the monomers. The neglect of the solvent does not affect 
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the static properties of chains in dilute solution. However, it does affect their physical 
dynamics (see Ref. 14). 

2 A Primer to Polymer Physics 
2.1 A Polymer in Good Solvent 

To substantiate the remarks of the introduction about the large-scale properties of polymers 
let us consider a specific example, a dilute solution of polyethylene. Polyethylene consists 
of CH2-monomers which are joined to form a Unear polymer (Fig. 1). A configuration of 
the chain may be specified by the positions of the monomers'* x = (ri, . . . , rV). Thermo- 
dynamic properties are calculated by averaging an observable A over all configurations 

(^) = 1 y Ax A{x) exp [ - /3C/(x)] . (1) 

Here (5 = k^T, Z is the partition function and U{x) the interaction potential. We assume 
that U{x) can be split into two parts:^ 

N-l 

U{x) = V C/o(&„..., 6,+^_) + C/i (a;, solvent) , (2) 
tt' ' ^ «' ' 

"short-range": i,6,ip,... "long-range" 

where bi = ri+i — ri denotes the bond vector from the ith to the {i + l)th monomer. 

The first term of Eq. (2), Uq, depends on the chemical nature of the polymer. It com- 
prises the potentials of the bond length £, the bond angle 9, the torsional angle (p, etc. 
(Fig. 1).'^ These potentials lead to correlations between the bond vectors bi and bj. Typi- 
cally, the correlations are of short range: they only extend up to some bond vector bj=i+i^ 

with in,a„ < N. 

Although distant monomers along the backbone of the chain are thus orientationally 
decorrelated, they can still come close in space. The resulting interaction is long-range 
along the chain backbone (Fig. 1). In Eq. (2), it is accounted for by the second term 
U\.^''' Ui depends strongly on the quaUty of the solvent.-^ In good solvents the monomers 
effectively repel one another (they want to be surrounded by solvent molecules), whereas 
they attract each other if the solvent cannot dissolve the polymer (bad solvent). 

Due to its long-range character, one intuitively expects C/i to influence the large-scale 
behavior of the chain more strongly than Uq. A possible test of this idea is to estimate how 
the size of a chain scales with TV. Common measures of the chain size are the mean-square 
end-to-end distance or the radius of gyration R"^ (Fig. 1) 

1 ^ 

R! = {{VN - nf) , ^g = ^ E ((^^^ - ^-)'> ' (3) 

1=1 

''Here, we adopt a description in terms of a so-called "united atom model". The united atom model repre- 
sents a CH2 -group by a single, spherical interaction site and does not distinguish between inner (CH2) and end 
monomers (CHs).'^ Furthermore, we neglect the momenta of the monomers to specify the configuration, as we 
assume the observables and interaction potentials to depend on positions only. 

'^Equation (1) does not contain the degrees of freedom of the solvent. They are assumed to be integrated out. 
Thus, U{x) is an effective potential -in fact, a free energy- depending on the properties of the solvent. 
■'^In Eq. (2) we assume that Uo is independent of the solvent quality. 
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Figure 1. Schematic illustration of polyethylene. The local properties of the polymer depend on its microscopic 
degrees of freedom: the bond length £, the bond angle 8, and the torsional angle 4>. Because the potential of 
the bond length is fairly "stiff", £ may be kept fixed at its equilibrium value £o in a modeling approach. By 
contrast, the potential of the torsional angle is much "softer". Thus, </», which characterizes rotations about 
a middle C-C bond, mainly determines the local conformation of the chain. All degrees of freedom (£,9,(p) 
determine the intrinsic stiffness of the chain. The stiffness reflects the persistence of orientational correlations 
along the backbone of the chain. Orientational correlations decouple on the length scale of the "persistence 
length" £p. For typical chain lengths, A'^ ~ 10^, £p is much smaller than the end-to-end distance iJc or the 
radius of gyration Rg. {Rg measures the average distance of a monomer from the center of mass (cm) of the 
chain.) Thus, the chain appears flexible on length scales larger than £p. If the polymer is dissolved in a good 
solvent, distant monomers (filled grey circles) repel each other when they come in contact. That is, the excluded- 
volume parameter u, measuring the effective interaction between distant monomers along the chain, is positive. 
Under these conditions (i.e., linear polymer with some flexibility and repulsive monomer-monomer interactions) a 
correspondence between the large-scale properties of the polymer and a critical system close to its phase transition 
can be estabhshed:*- may be identified with the reduced distance, t, to the critical temperature Tc of the 
phase transition, and or Rg scale with as the correlation length ^ of the order parameter does with r. u is 
a universal critical exponent, often called "Flory exponent" in polymer science. 



where R^m is the position of the chain's center of mass. Because i?e « Rg we focus on i?e 
in the sequel to illustrate the role played by Uq and Ui. 

Let bi denote the unit vector associated with the bond bi of fixed length £o (Fig. 1). 
Then, quite generally, we may write R^ as 



Apparently, the large-scale behavior of i?e depends on the range of orientational correla- 
tions between bond vectors. Two cases may be distinguished:^ 

®In part, the subsequent discussion closely follows that on p. 148 of Ref. 18. 



N-lN-1 N-1N-I~i 



(4) 
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1. If {hi ■ bi+k) is "short-range", i.e., if it decays more rapidly than 1/k for large k, the 
second term converges in the large- A?^ hmit. Then, 



= ml 



2^(6i-6i+fe)-l 



k=0 



=:Nel 



2^-1 



(N^oo), (5) 



where we introduce the persistence length £p in the last term, (^p measures the "persis- 
tence" of orientational correlations along the backbone and thus the intrinsic stiffness 
of the chain; see Fig. 1). Equation (5) shows that short-range orientational correlations 
only affect the prefactor -they renormalize the bond length tob — £o[2(^p/^o) ~ 1]^^^ 
(b is called "effective bond length"^)- but they do not change the scaling of i?e with 
N. The scaling is always "random-walk-like":'' i?e ^ N'^I'^J In polymer science, a 
chain exhibiting this random-walk-like behavior is commonly referred to as an "ideal 
chain". 

Of course, the finite-range correlations, assumed for in Eq. (2), are also of short 
range. Thus, provided Ui = 0, the end-to-end distance of a (long) chain is given 
by i?e = bN^^^, irrespective of the precise form of Uq. The microscopic degrees 
of freedom, i?, 6*, determine the prefactor, the effective bond length b, but not the 
scaUng with N. Therefore, if we are interested in studying large-scale properties, we 
can replace a chemically reaUstic model for polyethylene by a much simpler "coarse- 
grained model", which is microscopically unrealistic, but correctly reproduces the 
large- A'^ behavior. An example for such a coarse-grained model is a "bead-spring 
model", where N effective monomers ("beads") are coimected by harmonic springs 
of average length b (Fig. 2). 

However, if {bt ■ bi+k) decays as 1/k or more slowly (as l/k^ with y < 1) due to 
long-range correlations, the scaUng behavior of is changed. Instead of R^ N 
we find 

R', ~ ml fdk m . 6(0)) ~ I ^ j^'^ < Jj ; (6) 

Thus, long-range correlations lead to a "swelUng" of the chain size with respect to a 
pure random walk. 

Such long-range correlations are embodied in the potential Ui in Eq. (2). For a poly- 
mer in a good solvent a swelling of the chain dimension is physically reasonable. 
As soon as two (distant) monomers come close in space, they repel each other. On 
the level of the coarse-grained bead-spring model we can incorporate this repulsive 
interaction by writing Ui as (see e.g. Ref. 7 or the lucid discussion on pp. 16-20 of 
Ref. 16) 



N 

Ui{r^)= d?r-kBTvp{f)'^ with p(f') = V 5(r - n) . 



(7) 



''By the term "random-walk-like" we mean the diffusional motion of a Brownian particle. This motion can be 
thought of as resulting from the addition of many small displacements in random directions so that the overall 
mean-square displacement of the particle in time t, B?{t) = ([r(t) — r(0)]-^}, scales as ii^ ~ t. This allows for 
the following identifications in regard to polymer physics: i? <-> i?e and t <-> AT. 



5 




Figure 2. From a chemically realistic model to a coarse-grained bead-spring model. Local properties of the 
realistic model are determined by its microscopic degrees of freedom: 9, and </<. On the global level of the 
chain, however, the influence of the microscopic degrees of freedom can be lumped into one parameter, the 
effective bond length b. The microscopic degrees of freedom do not determine the scaling of the end-to-end 
distance, which, under the sole effect of Uq, is given by = bN^^^ ("ideal chain"). This behavior may 
be recovered from Eq. (1) when calculating Re with the potential Uq^ of a coarse-grained bead-spring model. 
This model identifies the monomers with spherical "beads" which are bovind to one another by harmonic springs 
with force constant Zk%T/h'^. (This bead-spring model is often called "Gaussian chain" model in the polymer 
UteratureJ) 

Here, p{r) is the monomer density at point r and v (> 0) is the excluded-volume 
parameter, v measures the strength of the repulsion of a binary contact between two 
beads. Because a binary contact occurs with probability p(r)^, Eq. (7) expresses the 
total energy penalty resulting from the repulsive contacts of all beads in the chain. 

From the previous discussion of Uq and Ui the following conclusion may be drawn: When 
focusing on the large-scale properties of Unear polymers with some flexibihty and pre- 
dominantly repulsive interactions we may forego a microscopic description in favor of a 
coarse-grained model. An example is the bead-spring model introduced above (Fig. 2), 
which is characterized by two parameters, h and v. Another possibility is a self-avoiding 
walk (SAW) on a (hyper-cubic) lattice. That is, a random walk which is not allowed to 
visit an already occupied lattice site again (see Sec. 4.1). The replacement "microscopic 
model SAW" is permissible because a linear polymer in good solvent can be shown to 
correspond to a critical system which undergoes a phase transition for A/^ ^ cxd (Fig. 1). It 
belongs to the universaUty class of the n- vector model in the limit n — > 0.^' This implies 
that the large- A?^ behavior is determined by critical exponents. For instance, 

i?g oc i?e = hN" or Z ~ N^-^ {N ^ oo) , (8) 

where the partition function Z counts the number of A^-step SAW's starting at the origin 
and ending anywhere. The connectivity constant /i and the bond length b are non-universal. 
They depend on the polymer and the external conditions (temperature, solvent, etc.). By 
contrast, the critical exponents i' and 7 are universal. They only depend on the dimension 
of space.* Thus, they can be determined for all polymers by studying this (simple) model. 

'In the coiffse of the research on critical phenomena it has become clear that all systems with short-range, isotropic 
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Figure 3. Schematic phase diagram of flexible polymers (see Chap. 9 of Ref. 16 or Chap. 4 of Ref. 23). For 
small monomer density p the solution is dilute. Three different regimes may be distinguished according to the 
temperature T: swollen chains [Eq. (8), T > Tq: dilute (I)], nearly ideal chains [Eq. (10), T Ri Tq: dilute 
(II)]. and collapsed chains [Eq. (11), T < Tq: dilute (III)]. There is an interval AT around the 0-point of 
order AT/Tq ^ l/\/~N, where the chains are nearly ideal. Whereas the chains may be considered as being 
isolated in dilute solution, they strongly overlap in the semidilute regimes. For T < Tc{N) phase separation in a 
dilute phase of collapsed chains and a semidilute solution of nearly ideal chains occurs. If the monomer density 
approaches 1, we obtain a polymer melt. At high T the melt is a (viscous) liquid, whereas at low T it may become 
a glassy^** or a semicrystalline-^ solid, depending on the ability of the polymer to form ordered structures or not. 



In fact, the currently most precise values of v and 7 (see footnote on page 30) have been 
obtained from high-precision Monte Carlo simulations of SAW's.^^-^^ 



2.2 Phase Diagram of a Polymer Solution 



The utility of coarse-grained models to investigate the statistical physics of polymer sys- 
tems is not limited to the previous example. A dilute solution in a good solvent is just 
one region in the phase diagram. The phase diagram of flexible polymers is schematically 
shown in Fig. 3. Out of the various regimes we choose to discuss two cases in more detail, 
a chain in another than good solvent and (high-temperature) polymer melts. In the follow- 
ing sections we concentrate on those cases because novel Monte Carlo approaches have 
been applied to them. 



interactions, the same dimension of space d, and the same dimensionality n of the order parameter (n = 1: scalar, 
n > 2: 71-dimensional vector) have critical exponents which depend only on (d, n) and take the same values as 
those of the n-vector model. 
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A Chain in a Q-Solvent or a Bad Solvent. To extend the discussion of the good solvent 
to other solvents let us reconsider Eq. (7). This equation corresponds to the first term of a 
virial expansion in the monomer density p{r) . That isj' 



Ui = / dV 



2 D 



(9) 



This identifies the excluded-volume parameter v with the second virial coefficient. In 
general, the virial coefficients depend on temperature T. The second virial coefficient 
vanishes at some temperature, called "6-temperature T@" in polymer science, and behaves 
as ?; = uo(l — Tq/T) close to the 9-point {vq = const. > 0). This impUes that we can 
tune the solvent quality by temperature. In addition to the case of a good solvent (T > Tq) 
two further cases may be distinguished: 

1 . Q-solvent (T = Tq): Since binary interactions are absent [but ternary interactions are 
present: w > in Eq. (9)], the polymer behaves nearly as an ideal chain: 

i?e oc i?g ~ \/N (+ In A'' corrections) . (10) 

2. Bad solvent (T < Tq): Since the binary interactions are attractive, the polymer is 
collapsed to a dense sphere of monomers, implying that the average monomer density 
p inside the sphere is of order 1. Thus, 

N_ . „ „ 1 

m 



=^ iJeOciJg^iV' with iy = - . (11) 



The simulation of this situation is complicated because the time to equilibrate the 
chain and to sample equilibrium properties from many independent configurations 
becomes exceedingly long. Two factors are responsible for that. On the one hand, 
the local dynamics is sluggish (maybe even glass-like) due to the dense packing of 
monomers that strongly attract each other. On the other hand, the polymer is entangled 
with itself. Bonds cannot pass through each other. These topological constraints may 
also lead to slow dynamics for long chains. 

The Size of a Chain in a Polymer Melt. In a good solvent a chain expands with respect to the 
ideal state, owing to long-range monomer-monomer repulsions. This is peculiar to dilute 
solutions. In a dense liquid of chains, a "polymer melt", the situation is quite different. 
One can show*'^'^^ that the intra-chain excluded-volume interactions are screened by the 
presence of the surrounding polymers. Thus, a chain in a melt behaves on large scales as 
an ideal chain, implying i?e oc i?g ~ N^^"^ (see Fig. 4). 

This ideaUty, first proposed by Flory,'^ appears fairly unexpected. Some feeling why 
this should be so may be obtained from the following argument: Inside the volume of a 
chain (~ i?g) the monomer density resulting from the N monomers of the chain is very 
small. For ideal chains it is of order N/R^ ~ N~^/'^, whereas it scales as ~ Ar-O-764 
under good solvent conditions (using Eq. (8) and jy = 0.588). We see that in dilute solution, 
swelling reduces the monomer density inside the chain and thus the total interaction energy 
[see Eq. (7)]. However, no energetic advantage may be gained in a melt because the overall 
monomer density is p ^ 1. Swelling would reduce the number of intra-chain contacts, but 
this reduction must be compensated by inter-chain contacts to keep p constant. Thus, a 
chain has to have N^^'^ contacts with other chains, which is huge in the large-iV Umit. 
This strong interpenetration of the chains suppresses the expansion of an individual chain. 
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Figure 4. End-to-end distance versus chain length A'^ for the (athermal) bond-fluctuation model which will 
be discussed in more detail in Sees. 4.1 and Sec. 5. Results for three volume fractions (of occupied lattice sites) 
are given, illustrating the dilute (<p = 0), the semidilute (0 = 0.03125) and the melt (<p = 0.5) limits of the 
schematic phase diagram (Fig. 3). Using the slithering-snake algorithm (Sec. 5.2) it is possible to simulate chains 
containing up to Af = 32768 monomers for ij> < 0.5. Since the slithering-snake algorithm becomes less efficient 
at high densities (Sec. 5.2), the recently proposed double-pivot algorithm, described in Sec. 5.3, was harnessed 
to probe systems of higher densities (4> = 0.5). Periodic boxes of linear size L = 512 and containing up to 
2-^^ monomers are required to eliminate finite-size effects. Such periodic boundary conditions are not needed 
for single chains. Here, an infinite box was used (L{<l) = 0) = oo). As only excluded-volume interactions 
are taken into account, good solvent statistics applies in dilute solution. The chains are swollen, as indicated by 
the exponent u = 0.588 (solid line), which fits the data over three orders of magnitude. In the opposite (so- 
called) melt limit long-range correlations appear to be screened down to small chain lengths of about N x 10 
(grey dashed line).^' Both chain statistics are visible for the intermediate density (c/) = 0.03125): Small chains 
{N <^ g, Rs <^ ^£ swollen (solid line) and long chains are Gaussian (dashed line). The intercept of both 
lines defines the size ^ of the "excluded volume blob"*-' and the number of monomers g that the blob contains. 
The indicated numbers are specific to the volume fraction (and persistence length) given, but are independent of 
chain length. For a given density § coiTesponds to the chain size where the coils start to overlap. Also presented 
in the figure is the spatial distance ((r; — along the longest chain for (f> = 0.03125 (dotted line). 

With the exception of small N ot k (i.e., N,k <^ 10) this distance is, within the numerical accuracy of the data, 
identical to iJc(Af) with N = k. This agreement also demonstrates that the difference between a segment of 
a long chain and a chain having the same length as the segment becomes irrelevant for distances larger than i;. 
In precisely this sense the (long-range) excluded volume interactions are screened in semidilute solutions and in 
melts. Mean-field descriptions become appropriate on the level of coarse-grained (Gaussian) chains of blobs. 

2.3 Dynamics of Polymer Melts: Rouse and Reptation Models 

The Rouse Model. As a monomer in a dilute solution moves, it creates a vortex ("wave") in 
the solvent. The solvent transports the "wave" which is transported to other monomers of 
the chains so that a coupling between the motion of (distant) monomers arises (see Ref. 14). 
This long-range hydrodynamic interaction becomes screened by other chains when the 
concentration of the solution increases.^ In a dense melt, hydrodynamic interactions are 
completely suppressed. Thus, it is generally believed that the Rouse theory^- ^ provides 
a viable description for the long-time behavior of polymer dynamics in a melt, provided 
entanglements with other chains, giving rise to reptation dynamics,^-^ are not important 



9 




chain topology 

Figure 5. Sketch of the reptation concept for the dynamics of long-chain polymer melts.^-' The chain is supposed 
to be enclosed in a "tube" formed by its neighbors. The tube may be chai'acterized by an axis, the primitive 
path. The tube confines the motion of the enclosed chain: It predominantly moves along the primitive path. 
Perpendicular excursions are suppressed beyond the tube diameter dj. The tube diameter is larger than the 
effective bond length b: dj = N^, where the "entanglement length" A^c 2> 1. The primitive path represents 
the shortest connection between the chain ends, which respects the topology imposed on the enclosed chain by the 
entanglements with its neighbors. The length L of the primitive path is thus larger than iJc, which is the shortest 
connection between the chain ends in space. L varies linearly with A^: Ldj = so that L = d-Y^N/Nc). 



(see Fig. 5 and also below).-' 

The Rouse theory assumes the chains to be ideal and models them as a sequence of 
Brownian beads, connected by harmonic springs and subjected to a local random force 
and a local friction.^- ^ This bead-spring model is characterized by two parameters: the 
effective bond length b and the monomer mobility m. The mobility, or more precisly 
1 /to, measures the time it takes a bead to diffuse over the distance b. Thus, the diffusion 
coefficient of a monomer is proportional to mb^. As the center of mass (CM) of a chain 
does not experience any external force other than the antagonistic friction and random 
forces, the theory predicts that the CM diffuses freely at all times 

53 (i) = {[Remit) - i?cm(0)]') = 6DNt , (12) 

where i?cm(i) denotes the position of the CM at time t. The diffusion coefficient of a chain 
is by a factor of N slower than that of a monomer, i.e., 

Dn^ — . (13) 

From Eqs. (12,13) the longest relaxation time tn can be obtained. Arguing that a chain is 
relaxed when its CM has diffused over a distance of the order of its own size, we find 

gsiTN) ~ Dntn Rl ^b'N ^ tm , (14) 



^To our knowledge, there is no established derivation of the Rouse model from a microscopic theory. For a recent 
attempt see Ref. 28. 



10 



where the ideality of the chain was exploited. 



Strongly Entangled Polymers and Reptation Model. The single-chain picture proposed by 
the Rouse theory is supposed to be vaUd as long as entanglements with other chains do not 
dominate the polymer dynamics. This is beUeved to be the case for short chains, for which 
N is smaller than the entanglement length (Fig. 5). For S> the prevailing picture 
is that a chain is enclosed in a temporary "tube" formed by its neighbors. Entanglements 
force the enclosed chain to diffuse along the contour of the tube having a length of L N 
("reptation"; see Fig. 5).^'^ Because the curvilinear diffusion through the tube is presumed 
to be Rouse-like, reptation theory predicts the relaxation time of the chain to scale with A'' 
as 

so that the diffusion coefficient of the CM in space is given by [Eq. (14)] 

93{tn) ~ Dntn ^Rl^N Dn^^ . (16) 

Experimentally, one finds a stiU stronger dependence: tjv ^ N^^-"^ and ^ ]y~-2.4k 
Clearly, simulation methods which attempt to model the true physical dynamics, such 
as Molecular Dynamics or Monte Carlo algorithms employing local random moves, must 
suffer from these long relaxation times. Various alternative Monte Carlo methods have 
been proposed to efficiently equilibrate dense polymer melts. We will present some of 
these approaches (Sec. 5 and Sec. 6). 



3 Monte Carlo Methods: A Brief Overview 

In equilibrium statistical mechanics thermodynamic properties are calculated as ensemble 
averages over all points a; in a high-dimensional configuration space F.' In the canonical 
ensemble the average of an observable A{x) is given by 

{A)= j 6x A{x)P^{x) = ^ J dx A{x) exp [ - l3U{x)] . (17) 

In general, the integral cannot be solved analytically. Monte Carlo (MC) simulations pro- 
vide a numerical approach to this problem by generating a random sample of configuration- 
space points xi, . . . , Xm, ■ ■ ■ , xm according to some distribution Ps{x). (A) is then esti- 
mated by^^"^^ 

M M 

= . (18) 

m=l m=l 



*^This exponent varies very little -if at all- with the chemical properties of the (Unear) polymer.^'-''* 

'We assume that the momenta can be integrated out, since the observables only depend on the positions of the 

particles. 
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where we introduced the "weight" W{x) = Peq{x) /Ps{x). Note that, while (A) is a num- 
ber, A is still a random variable. Whether A represents a good estimate for {A) depends on 
on the total number M of configurations used and, for a given M, on the choice of Ps{x). 
To see this in more detail™ let us define the mean value with respect to Ps by 



((•))3 = Jdxmix). 



(19) 



34 



For large M, the average of A and its variance vais{A) may be estimated from the small- 
fluctuation approximations 

'Y 



vars 



This gives" 
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{Y)l{Z'). 
{Z)l 
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1 

M 



W{A-{A)y 



(20) 



(21) 



(22) 



Equation (22) shows that A provides an unbiased estimate of {A) in the limit Af 1 
unless {W) ^ 1, i.e., unless Ps(ic) is very different from P^n{x). When Pf,{x) deviates 
considerably from Peq(a^), it predominantly samples configuration-space points, which are 
not representative of the thermal equilibrium. One could try to compensate this inefficient 
sampling by making M larger and larger. However, on the one hand this renders the simu- 
lation very time-consuming. On the other hand, there is no guarantee that the maximum M 
one is willing (or able) to simulate suffices to outweigh the error incurred by the inadequate 
choice of Ps{x). 

Thus, Ps{x) should approximate Ptq{x) as closely as possible to obtain meaningful 
results from MC simulations. To this end, two approaches may be distinguished:^^- 

1 . Static MC methods: Static methods generate a sequence of statistically independent 
configuration-space points from the distribution Ps{x). In this case one has to tune 
the algorithm cleverly so that the weights W{x) do not get out of hand. Examples 
how to achieve this will be discussed in Sec. 6. 

2. Dynamic MC methods: Dynamic methods generate a sequence of correlated con- 
figuration-space points via some stochastic process which has P(,q{x) as its unique 
equilibrium distribution. In practice, this process is always taken to be a Markov 
process. ^^'-^^ The defining property of a Markov process is that it has no "memory". 
That is, the probability for the occurrence of the future configuration x depends only 
on the present configuration x' and not on the other configurations that the process 
visited in the past. 

Dynamic MC methods have become a widely used simulation technique, to which we will 
also heavily refer in the foUowing sections. So, we provide a brief introduction here (many 
more details may be found in Ref. 36). 



'"In part, our discussion closely follows Sec. 2.3 of Ref. 34. 

"Note that (W)s = 1, {WA)s = (A), etc. In Eq. (22) the Ri-sign means that there are corrections of ©(l/M^) 
which we have neglected. 
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Let us assume that the configuration space is discrete and that the Markov process 
evolves in this space in discrete time steps At (= 1). The time evolution of this Markov 
chain may be characterized by the "master equation" for the probability P{x, t) to find the 
system in the state x at time t 



P{x,t+1) - P{x,t) = \^w{x\x')P{x',t) -w{x'\x)P{x,t) 



(23) 



Here, w{x\x') denotes the transition probability from x' to x which is independent of 
time. (In the continuous time limit (At 0) it becomes a "transition rate", i.e., a transition 
probabiUty per unit time.) Equation (23) expresses the balance between the flux of all other 
states x' towards x (first term of the rhs), leading to an increase of P{x), and the flux away 
from X (second term of the rhs) which diminishes P{x). Note that only terms with x ^ x' 
contribute. We can rewrite Eq. (23) by including the missing term for x = a;' if we 
take into account that w{x\x'') is normalized. Since a transition from x' to some state x, 
including x' , will occur with certainty, w{x\x'^ satisfies 

Yw{x\x') = \ . (24) 

a:' 

Inserting Eq. (24) into Eq. (23) the master equation takes the following form 

P{x,t^\)=Yw{x\x')P{x' ,t) . (25) 

a;' 

For the application of these results to statistical physics it is necessary that -P(a;, €) con- 
verges to a unique stationary distribution, irrespective of the initial configuration of the 
system, in the long-time limit and that this distribution is the (canonical) equilibrium dis- 
tribution (x) . Thus, the right-hand side of Eq. (23) must vanish for P{x' , t) = P^ {x'). 
Certainly, this is the case if each term of the sum vanishes separately. This leads to the con- 
dition of "detailed balance" (see Refs. 31-33, 36) 

w{x\x')P,^{x') = w{x'\x)P,^{x) . (26) 

To exploit this condition in MC algorithms the transition probability may be split into 
two independent parts: First, we propose a transition from x' to x according to some 
probability Ppro{x' x). Then, this move will be accepted or rejected with probabilities 
acc(a;' x) and 1 — acc(a;' — > x), respectively. So, we have 



w{x\x') ^ Pp,o{x' x)acc{x' x) ^ 
w{x'\x) Ppro{x — > x') acc(x — > x') 

To solve this equation for acc(a;' — > x) we set 

Pp,„(x^x')e-''^(«') 



(27) 



.cc(x-^x) = f^ ^-,^^;^_^„„,, J. (28) 

From Eq. (27) we see that the function F{x) satisfies F{x) /F{l/x) = x. One solution to 
this equation was proposed by Metropolis et al.-?^ F{x) = min(l, x). This leads to the 
"Metropolis criterion" for the acceptance probability 

acc{x' ^x)= min ( 1, ^P^('^^ ^-f3[u(.)-ui.')]\ (29) 
V ^pro(a:;' ^ x) ) 
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The Metropolis criterion is the core of essentially all dynamic MC algorithms. It embodies 
detailed balance which guarantees that the simulation, irrespective of the initial configura- 
tion, converges to the canonical equilibrium distribution, provided the set of chosen Monte 
Carlo moves leads to ergodic sampling." 



Detailed Balance versus Stationarity. Detailed balance is an important, but very strict 
criterion. Less stringent is the condition of stationarity [Eq. (25)] 

) = ^«;(x|a;')Peq(a;') , (30) 

a; 

implying that Peq{x) remains invariant under the Markov dynamics. Stationarity in con- 
junction with the ergodicity of chosen set of MC moves ensures a valid simulation. 

In practice, this milder condition may be important. Imagine that we want to update 
a polymer chain consisting of N monomers and that each monomer can be displaced in 
iVdis directions. One possibility is to select a monomer and a direction randomly. Thus, 
Ppro{x' x) = l/{NNais) = PpTo{x — > x'). This procedure obeys detailed balance: In 
the next move the same monomer and the reverse displacement may be chosen with the 
same a priori probability. On the other hand, one could also attempt to move one monomer 
after the other, proceeding regularly from monomer 1 to monomer N. This sequential 
updating scheme violates detailed balance: The next step never selects again the monomer 
whose displacement has just been attempted. So, the probability for the reverse move is 
zero. 

However, sequential updating is a valid scheme if the individual steps obey detailed 
balance"*^ or at least stationarity. To see that we can write the transition probability from 
x' to X for sequential updating as 

M;(a;|a;') =^---^^w(^Ha;|zAr)---w;(2)(z2|2i)u;('H2i|a;') ■ (31) 

Zn Z2 Zl 

This means that the process passes sequentially first with probabiUty w^^^ {zi\x') from x' 
to Zl by attempting to move the first monomer, then from Zi to 2:2 by attempting to move 
the second monomer, and so on until configuration x is reached. Multiplying Eq. (31) by 
Peq{x') and summing over all x' we find 

Y,wix\x')P,^{x') 

x' 

= E • • • E E ^^'^^('^l^iv) • • • W^'\Z2\Z^) W^'\Z,\X')P,^{X') 

Zn Z2 Zl ^ V ' 

= feq(2l) 

= . . . = Peq(a;) , (32) 



°By "ergodic" sampling we mean that the probability of finding the system in configuration x, starting from some 
state as' (including as), is non-zero for all x after a sufficiently long time.^* This definition is a bit dangerous 
because it conflicts with others in the Uterature. For instance, in mathematical texts on Markov chains (= discrete- 
time Markov processes with a discrete configuration space) our definition rather corresponds to an "irreducible 
and aperiodic" chain (there, "ergodic" is a synonym for "irreducible").'''' In Ref. 40 our definition would be 
termed "regular sampling". 
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i.e., sequential sampling preserves the stationarity of the equilibrium distribution. Thus, 
it is a correct simulation procedure. This conclusion is important for a variety of MC 
methods which perform different trial moves in a fixed order. 

4 Some Coarse-Grained Simulation Models 

In Sec. 1 we introduced the term "coarse-grained model". This was defined as a model 
which associates a group of chemical monomers with a "bead" (effective monomer) in 
order to eliminate microscopic degrees of freedom (bond length vibrations, etc.). Here, we 
refine our definition and distinguish between two types of coarse-grained models: 

1. The coarse-grained model is derived from a specific polymer In practice, this usually 
implies that the properties of the model (potential parameters, density, etc.) have to be 
adjusted to results from atomistic simulations of the polymer under consideration (see 
Appendix A for an example). The incentive to devise such models rests upon the fact 
that they may be simulated much more efficiently than their atomistic counterpart. 
Thus, it is tempting to split the simulation into two levels: First, one uses the coarse- 
grained model for equiUbration and for the determination of large-scale properties. 
Then, atomistic details may be reinserted to allow for a thorough comparison with 
experiments. Recent attempts to perform such multi- scale approaches are described 
in Refs. 41, 42 (see also Appendix A). 

2. The coarse-grained model has no direct connection to any specific polymer It is 
a generic model retaining only features common to all polymers of the same chain 
topology. For (uncharged) linear polymers these features are chain connectivity, 
excluded-volume interactions, and, additionally, monomer-monomer attractions if 
one wants to simulate O- or bad-solvent conditions (see Fig. 3) . Many of these generic 
models, be it lattice or continuum models, have been introduced in the literature (see 
Refs. 43, 44 for comprehensive overviews). In the following we present those models 
in more detail, which will be discussed in Sees. 5,6. 

4.1 Lattice Models 

The Self-Avoiding Walk. About 50 years ago Orr and Montroll"** proposed the self-avoiding 
walk (SAW) as a model for a linear polymer in a good solvent. The SAW is defined 
on a discrete lattice, often on a square or simple cubic lattice (Fig. 6). Each monomer 
occupies one lattice site, the bond length equals the lattice constant, and the bond angles 
are restricted by the lattice geometry and by the repulsive hard-core monomer-monomer 
interaction (e.g. 90° and 180° for the cubic lattice, as immediate backfolding is forbidden). 
This model can be complemented by attractive interactions if, for instance, an energy gain 
— e is associated with every occupied nearest neighbor pair.'*^ In addition to excluded 
volume interactions the simulation then also has to take account of the Boltzmann factor 
exp(nnne/fcB^r), where rinn is the number of nearest neighbors. 

To simulate the SAW by dynamic Monte Carlo one must first decide about the ele- 
mentary moves that propose a new SAW configuration x from an old one x'. The earli- 
est suggestion"'^ comprised one-bead excitations^'*'''^'^ (Fig. 6). In these algorithms, one 
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kink jump 



end rotation 



crankshaft 



jump forbidden 



Figure 6. Left figure: Single-site self-avoiding walk (SAW) of chain length TV = 10 on a simple cubic lattice 
(solid lines and black dots). The grey dots and the grey dashed lines indicate the moves discussed in the text: 
end-bond rotation, kink jumps and 90° crankshaft motion. Right figure: Sketch of a possible configuration of 
monomers in the 3D bond-fluctuation model (BFM). (A vectorized version of the BFM algorithm can be found in 
Ref. 45.) The bond vector (3, 0, 0) (thick black aiTow) blocks four lattice sites (marked by Q) that are no longer 
available to other monomers due to the excluded volume interaction. This interaction also prevents the jump of 

the grey monomer in the direction of the large arrow ( >), since the comers of the monomers, indicated by 

would then occupy the same lattice site. 



chooses a monomer at random. If the monomer is at the chain end, the bond to its neighbor 
is turned to a randomly selected lattice direction. Due to the fixed bond length an inner 
monomer is only mobile if its bond angle is 90° on the square or simple cubic lattice. In 
this case, one attempts a "kink-jump" motion, i.e., a one-bead flip to the opposite lattice 
site. End-bond rotation and kink jumps are accepted according to the Metropolis criterion 
if the target sites are empty. 

These moves are special examples of the class of "local A^-conserving moves". Quite 
generally, a "local move" alters the configuration of a small piece of the original SAW 
while leaving the remaining monomers unchanged. This definition opens the possibility to 
invent moves comprising more than one bead, such as two-bead or three-bead excitations. 
Figure 6 shows a common example, the 90° crankshaft motion (only possible in 3D). The 
crankshaft motion removes an important drawback of kink jumps. It introduces new bond 
vectors, whereas a kink jump does not. Therefore, if only end-bond rotations and kink 
jumps are allowed, new bond orientations have to diffuse from the ends toward the interior 
part of the chains. This algorithm is not very efficient in reshuffling the bond vectors and 
so in preparing independent configurations. The inclusion of crankshaft motions remedies 
this problem. 

However, even then a disturbing feature remains. It has been proved that all local A^- 
conserving algorithms for two- and three-dimensional SAW's are not ergodic for large TV."*^ 
There are dense configurations ("double cul-de-sac" in 2D, "knots" in 3D; see Ref. 34) 
which are completely frozen: They can neither be transformed into nor reached from other 
configurations. Whether this problem is serious in practice is a question that, to our knowl- 
edge, is not fully settled (see e.g. Ref. 43 or footnote 9 of Ref. 50). One can argue that, 
if one starts from an extended configuration -for instance, from a straight rod- and if one 
is interested in high-T properties only, non-ergodicity effects due to compact structures 
should be small. This argument may be true for short chains,^' but should fail for long 

J'Here, it is not clear what "shorf really means. For N < 10^ the error incurred by using local A'^-conserving 
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ones, since it has been proved that the fraction of SAW's belonging to the ergodicity class' 
of the straight rod is exponentially small in the large- A?^ limit. Of course, if one is interested 
in low-T properties, problems with non-ergodicity might be sizable, even for small chain 
length.^' 

The Bond-Fluctuation Model. The bond-fluctuation model (BFM) was proposed^^'^-' as an 
alternative to a (single-site) SAW model, which retains the computational efficiency of the 
lattice without being plagued by severe ergodicity problems. The key idea is to increase 
the size of a monomer which now occupies, instead of a single site, a whole unit cell of the 
lattice (e.g. a square for the 2D- or a cube for the 3D hyper-cubic lattice; see Fig. 6). This 
enlarged monomer size has two important consequences: 

1. A priori, many different bond vectors can occur. This multitude is restricted by two 
conditions. First, adjacent monomers may not overlap. This limits the bond length 
to ^ > ^min = 2 (in units of the lattice constant). Second, the hard-core monomer- 
monomer interaction should suffice to prevent two bonds from intersecting each other 
in the course of the simulation. In 2D this only imposes an upper bound on the bond 
length, £ < 

^max — v^TS,^^'^^ whereas in 3D, in addition to £ < ^max — V'lO, some 
smaller bond vectors also have to be excluded.^^ The resulting sets of allowed bond 
vectors are: 



where [ • ] denotes a class of bond vectors sharing the same length, but differing in 
direction. For instance, the class [2, 0] ([2, 0, 0]) comprises all vectors with a length of 
2 and direction along the lattice axis (4 directions in 2D, 6 in 3D). Equation (33) gives 
rise to 41 bond angles in 2D^^ and to 87 bond angles in 3D.^'* This has to be compared 
to 3 (2D) or 5 (3D) bond angles for the SAW model on the hypercubic lattice where 
a monomer is associated with a lattice site. Due to the multitude of different bond 
lengths and bond angles the BFM is much closer to continuous-space behavior than 
the single-site lattice modeF.^^ 

2. Ergodicity problems are much less severe than for the single-site SAW. For the BFM a 
local iV-conserving move consists of selecting a monomer at random and of attempt- 
ing a displacement by one lattice constant in a randomly chosen lattice direction. 
As these local jumps* permit transitions between different vectors, the algorithm can 
escape from configurations where a single-site model would be frozen in.^^ If the 
attempted displacement satisfies both the bond vectors constraints of Eq. (33) and 
the excluded volume interaction, the move is accepted. Of course, it is also possible 

algorithms seems to be small (see Ref. 43 and the footnote 9 of Ref. 50). 

'By "ergodicity class of a straight rod" we mean all mutually accessible configurations, one of which is the rod. 
''The main advantage of lattice models is their computational efficiency. Longer length and times scales may be 
probed. However later on, the results of the simulation shall be compared to theories or experiences, which "Uve" 
in continuous space. So, the important question arises of how well the lattice algorithm approximates continuum 
properties. A general, intuitive answer is: The finer the lattice, i.e, the more sites are occupied by one particle, the 
closer the continuum limit should be realized. Recently, this statement was made more precise by the example of 
monatomic fluids interacting via a Lennard- Jones or a Buckingham potential.'* 
"Larger jvmips distances were also tested (in 3D), but fovmd less efficient in concentrated solutions.''* 



{&} 



{6} = [2, 0] , [2, 1] , [2, 2] , [3, 0] , [3, 1] , [3, 2] (2D) , 

[2, 0, 0], [2, 1, 0], [2, 1, 1], [2, 2, 1], [3, 0, 0], [3, 1, 0] (3D) , 



(33) 
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to include a finite interaction energy. Then, the move is accepted according to the 
Metropolis criterion. A possible choice* is to work with an energy — e between pairs 
of monomers with distance 2 < r < \/6- This interval comprises all neighbors which 
contribute to the first peak of the pair-distribution function^^ in a dense polymer sys- 
tem. This choice was made in studies of the G-point*" and of the phase transition in 
binary polymer blends (see Ref. 11). 



4.2 Continuum Models 

Two Bead-Spring Models. A widely used continuum model is the bead-spring model in- 
troduced by Grest and Kremer.^^ In this model nearest-neighbor monomers along the 
backbone of the chain are bonded to each other by a FENE (finitely extendible non-linear 
elastic) potential 



f 2 ^-^max [l (^/-^max)^] ^ < -^max j j-^^-j 

[ oo else , 



whereas all monomers, bonded and non-bonded ones, interact via a truncated and shifted 
Lennard-Jones (LJ) potential 

f/ts ^ r 4e [(a/r)i2 - [a/rf] + C{r,^,) for r < r^ut , (35^ 
1 else ^ 

where C(rcut) ensures that the potential vanishes at the cut-off parameter rcut- Such a cut- 
off is commonly employed to render the interaction short-ranged (Fig. 7)." The parameter 
e defines the energy scale and a the length scale of the system. That is, we set e = ct = 1 
(LJ units) in the following. 

For small values of the bond length the FENE potential is harmonic ("elastic behav- 
ior"), i.e., Up{i) = for < £ <C £max> whereas the logarithmic divergence imposes 
^ < ^max ("finite extensibility"). The parameters Iroax and k have to be chosen such that 
the possibility of bond crossing becomes so unlikely that it never occurs. Reference 61 
suggests fc = 30 and ^max = 1.5 (in LJ units). This has become a standard choice. 

The FENE potential alone does not prevent monomers from overlapping. To real- 
ize excluded volume the LJ-interaction has to be taken into account also between bonded 
monomers. The superposition of the FENE- and the LJ-potentials yields a steep effective 
bond potential with a minimum at £0 — 0.96 (Fig. 7). The shape of the bond potential 
depends on the cut-off parameter of the LJ-interaction: 

• If one takes rcut = ^min = 2^/^ (C(T"cut) = 1). i-c, as the minimum of the LJ- 
potential, the monomer-monomer interaction becomes purely repulsive. This model 
is commonly called "Kremer-Grest model". For isolated chains it realizes good 
solvent conditions. 



*Another choice uses a discretization of the Leimard- Jones potential.^^ 

"From a computational point of view short-range interactions are convenient because the simulation can be 
speeded up by neighbor lists.^'-*^ However, as the truncation ignores the contribution of the tail of the potential, 
the error incurred must be corrected before comparing with results for the full potential. For instance, the trunca- 
tion shifts the location of the critical point of the Uquid-gas transition in a LJ-Uquid (see Ref. 63 or Sec. 3.2.2 of 
Ref. 31 for details). To avoid these truncation effects some authors prefer to work with the fuU LJ-potential.*'* 
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Figure 7. Bond and Lennard-Jones potentials versus the distance r between two monomers of the bead-spring 
model (for the bond potential r = b). The bond potential results from the superposition of Eqs. (34,35). For both 
cut-off parameters rcut the bond potential was shifted by —20 to show it on the same scale as the LJ-potentials. 
The LJ-potential with rcut = 2^/^ is purely repulsive, whereas the potential with rcut = 2x2^/^ has an attractive 
minimum at rmia = 2^/^. 



• The simulation of 0- or bad solvents requires to incorporate part of the LJ-attraction 
by increasing rcut- Obviously, there is freedom where to cut off the attractive part. One 
possibiUty is rcut = 2 x 2^/^ (C'(rcut) = 127/4096).*^ This choice is a compromise 
between the wish to include the major part of the attractive interaction and the need to 
keep the potential short-ranged. The resulting phase diagram was studied in Ref. 66. 

(Yet) Another Bead-Spring Model. If we recall the idea of the coarse-graining -a coarse- 
grained monomer stands for a group of chemical monomers- it appears plausible that 
coarse-grained monomers are softer than their chemical counterparts. Thus, an exponent 
smaller than 12 in Eq. (35) may be better suited to represent their repulsion. In fact, such an 
observation was made in a recent effort to develop a coarse-grained model for poly(vinyl 
alcohol) (see Ref. 67 and Appendix A). This study also suggests the following generic 
model which may be considered as a variant of the Kremer-Grest model. 

In this ("Kremer-Grest-Uke") model non-bonded monomers interact via a purely repul- 
sive 9-6 LJ-potential, 

C/'^P(r) = ( ^0 [(^o/r)9 - {cJo/rf] + C{rr^) for r < r„^„ = (3/2) V3 , ^^^^ 
I else ^ 

where eo = 1-511 and C{rmin) = 4eo/27. These non-bonded interactions are excluded 
between nearest neighbors in the chain, which are connected to each other by a harmonic 
potential 



(A; = 2141.84 C7^^ 4 = 0.97ao). 



(37) 



The equiUbrium bond length £q agrees with that of the Kremer-Grest model. The spring 
constant k has to be chosen so large to inhibit bond crossings (see Ref. 68 for further 



19 



discussion). A similar bond potential, in conjunction with Eq. (35) and rcut = 2^/^, has 
recently been used to study the effect of the bond length on the structure and dynamics of 
polymer melts.^' 

Local Moves for Continuum Models. The continuum models are constructed for use in 
Molecular Dynamics simulations. However, simulation within Monte Carlo schemes is 
also possible. Similarly to the lattice models a local updating scheme can be realized by 
selecting a monomer and a direction at random and by attempting a displacement in the 
chosen direction. This proposition is again accepted according to the Metropolis criterion. 

The size A of the displacement is a tunable parameter. It should neither be too small 
nor too large. If A is too small, many moves may be accepted, but the system advances 
only slowly in configuration space. Many displacements are thus needed to obtain well 
decorrelated configurations. On the other hand, if A is too large, many moves will be 
rejected and the decorrelation is also slow. A scheme how to optimize the choice for A is 
explained in Sec. 3.3 of Ref. 31. 

5 Monte Carlo Methods for Polymers: From Local to Non-Local 

Moves 

The method of importance sampling is based on a Markov process in configuration space. 
A priori, this stochastic dynamics is merely a numerical algorithm, aiming at an efficient 
sampling according to Peq{x). It need not correspond to the physical dynamics of the 
(polymer) system under consideration. An appealing consequence of this feature is the 
freedom to invent clever MC moves which decorrelate the configurations in the smallest 
(CPU) time possible. These non-physical moves serve to rapidly equilibrate the system 
and to produce statistically independent equilibrium configurations for the study of struc- 
tural and thermodynamic properties. We will pursue this idea in Sees. 5.2,5.3. In the 
following section we rather want to concentrate on local moves and the ensuing dynamic 
interpretation of the MC method. 

5.1 Local Moves: Studying Dynamic Properties with Monte Carlo 

By employing non-local moves we can explore the statics of the system, but information 
about its dynamic properties is lost. Of course, the equilibrated configurations could be 
used in a Molecular Dynamics (MD) simulation to analyze the dynamic properties. How- 
ever, if one is not willing to do that, the question arises of under which conditions the MC 
dynamics can be reahstic. The answer to this comprises two parts: 

1. Certainly, one can only expect the MC dynamics to become reUable on length and 
time scales where the deterministic motion of the monomers has been damped by the 
interaction with the surrounding (other monomers and/or solvent). For instance, in a 
(classical) MD simulation the monomers move ballistically at early times, i.e., their 
displacement is proportional to t. This is a consequence of the underlying Newtonian 
dynamics in the limit of vanishing force. At short times the monomers behave as 
if they did not "feel" the bonding to their neighbors and the presence of other par- 
ticles, that is, as if they were free particles. As time increases, the interaction with 
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the surrounding becomes important. The motion of the monomers is then a result 
of a multitude of individual collisions. This "averaging" over fast degrees of free- 
dom gradually lends a stochastic character to the dynamics which ultimately becomes 
diffusive in the long-time limit. 

2. The moves should be "physical". Usually, this implies that they are local. ^ Further- 
more, the dynamics should not be dominated by the momenta which are absent in 
Monte Carlo. The latter condition is satisfied in dense melts, but not in dilute solu- 
tion. In dilute solution the motion of distant monomers along the chain backbone are 
coupled via hydrodynamic interactions (see Ref. 14). and Sec. 2.3). Thus, we might 
expect that a local Monte Carlo algorithm reproduces Rouse dynamics where these 
long-range interactions are neglected. 

This expectation can be verified by estimating the scaling of local MC algorithms with 
N. To this end, let us assume that the center of mass (CM) of an isolated chain, be it 
on a lattice or in the continuum, may be considered as a free Brownian particle. This is 
reasonable, since the CM does not experience any external force other than the random 
force of the heat bath (resulting from the compound effect of the random monomer hops 
and the acceptance criterion). So, it should diffuse freely [Eq. (12)]. The corresponding 
diffusion constant Djv depends on chain length. To estimate this dependence we can argue 
that the center of mass is displaced by ~ 6/-/V, if one monomer moves over a distance of 
order b while the other monomers remain fixed. This elementary motion takes on average 
the time 1 /m with m denoting again the (temperature, density, etc. dependent) mobility 
of the monomer. For the CM to diffuse over the distance b, N such random motions are 
needed. This take the time mx N, which we use as our time unit here.'" Utilizing Eq. (12) 
we then find 93(1 = 1) ~ (mN) x {b/Nf ~ Dn. So, 

Inserting this result in Eq. (14) we obtain the relaxation time of a chain 

^'^^+^'^ r iV*2 (jjjgjji chain: v = 0.5) , 
~ TO " \ TV^s.ire (3j) excluded-volume chain: v = 0.588) . ^ ' 

Equations (38,39) agree with the predictions of the Rouse theory [Eqs. (13,14)]. 



Monte Carlo Dynamics versus Molecular Dynamics: An Example. The previous argu- 
ments suggest that the MC dynamics, based on local moves, becomes reaUstic for time and 
length scales outside the microscopic regime (of a bond). We want to support this assertion 
by a comparison between MC and MD simulations. 

Figure 8 shows the diffusion coefficient Dat of a chain versus chain length. Dn was 
derived from the long-time limit of Eq. (12) for both the BFM, simulated via MC, and 
the Kremer-Grest-like continuum model of Eqs. (36,37), simulated via MD. The figure 

"Examples for local moves of lattice models are given in Fig. 6. See the very end of Sec. 4.2 for a brief discussion 
of local moves in continuvmi models. 

'"This statement introduces the time unit Tmcs of a "Monte Carlo step (MCS)". A MCS is defined as the time 
it takes to give each of the N monomers the possibiUty to move once.'^''^ Thus, we measure time in units of 
attempted elementary moves per monomer. 
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Figure 8. Diffusion coefficient Df^ versus A'. Two simulation metfiods are compared: Tlie open symbols rep- 
resent MC data of the (lattice) BFM, the filled symbols were obtained from MD simulations of a (continuum) 
Kremer-Grest-like model [Eqs. (36,37)]- For both models results obtained in 3D for a dilute solution and for a 
melt are shown. For the BFM this corresponds to the following volume fraction <j>: <p = 0.0078 (dilute), (j> = 0.5 
(melt). For the Kremer-Grest-like model this corresponds to the following monomer densities p: p = 0.0835 
(dilute), p = 0.835 (melt). Qualitatively, the MD simulations yield the same dependence of Djv on A'. To 
illustrate this agreement the MD data were vertically shifted by an amount that optimizes the agreement with the 
MC results. (The shift factors are different for the dilute solution and the melt.) In dilute solution, we find Rouse 
behavior [Eq. (13)] for both methods. In the melt, the chains diffuse more slowly. The dependence of Djv on 
N is qualitatively compatible with the Rouse-to-reptation crossover when N passes the threshold (Sec. 2.3). 
Quantitatively however, there are deviations. Particularly for large N, the decrease of Dm appears to be stronger 
than predicted by reptation theory [Eq. (16)]. Roughly, we find ~ AT**"^-*. Adapted from Ref. 70. 



displays the results of the simulations for a dilute solution and a dense melt.^ Clearly, 
there is a high degree of accord between the results, illustrating that the BFIVI with local 
moves reproduces the reahstic dynamics of the IVID simulations. Thus, IVIC simulations 
can be more than just a versatile tool to calculate high-dimensional integrals. They may 
provide information on the dynamics of a system.^ 



^In the BFM, density is commonly specified in terms of the volume fraction (p of lattice sites occupied by 
monomers. As a monomer comprises all sites of a unit cube, the monomer density p is smaller than by a factor 
8, p = 0/8. Although the value (p = 0.5 appears small, the work by Paul et al?'' established that the chains 
have melt-like properties at this density (see also Ref. 4). Since then, (p = 0.5 has become a standard choice (in 
3D). For the Kremer-Grest model, the work of Ref. 61 showed that a monomer density of p = 0.85, or a value 
close to this, is a good choice to realize melt-Uke behavior. We adopted this choice in our MD study. The MD 
simulations were done at constant volume and constant temperature (Langevin thermostat*^). 
^There is further ample evidence for the correctness of this statement from other studies. For simple liquids of 
LJ-particles see e.g. Ref. 71. For polymers see the review in Ref. 8 or the comparison of MD simulations for 
polybutadiene and polyethylene with MC simulations of the BFM.^^ Furthermore, Monte Carlo methods have 
been applied to simulate dynamic processes in such diverse fields as relaxation phenomena in spin and structural 
glasses, spinodal decomposition of mixtures, nucleation processes, diffusion-limited aggregation, etc. (see e.g. 
the textbooks of Binder and Heermann'^ or Landau and Binder'^). 
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Figure 9. Slithering-snake move (a) and general reptation moves (b,c). Both moves are illustrated by the 
shrinkage-growth implementation. For the slithering-snake algorithm, a randomly chosen end bond (dashed 
line) is removed and then a new bond vector (also randomly chosen) is attached to the other chain end. For the 
general reptation algorithm, thi'ee moves are shown: Kink-kink transport (b) and kind-end/end-kink reptation (c). 
Kink-kink transport implies that a randomly chosen kink is shrunk to a bond and a new kink is inserted some- 
where else along the chain. Kink-end reptation (—*) amounts to replacing a randomly chosen kink by a bond 
and to appending two new bond vectors (also randomly chosen) to the other chain end. End-kink reptation (<— ) 
corresponds to the reverse "reaction". 



Relaxation Time and Computational Complexity. An important issue in any algorithm is 
its "computational complexity". Quite generally, the computational complexity may be 
defined as the time required to solve a computational problem.^'' Here, the computational 
problem is to decorrelate chain configurations. According to Eq. (39) this takes a relaxation 
time tm ^ N^^"^" in units of the Monte Carlo step (MCS; see footnote on page 21). As 
a MCS comprises N attempted moves of a monomer, the computational complexity Tcc 
scales with N as Tcc = Ntn ^ N^+^''. 

This rapid increase of Tcc with chain length -called "critical slowing-down"-'-'- makes it 
difficult in practice to efficiently decorrelate configurations of long chains by local moves. 
In order to be able to simulate large chains with sufficient statistics, moves have to be 
implemented, which reduces (t^c ~ N" with a < 2 + 2v) or even eliminate (Tcc ~ iV°) 
the critical slowing-down. These moves cannot be local, they have to act, in some way, 
on all monomers of the chain. In the following we want to discuss two examples of such 
global updates: bilocal moves and the pivot algorithm. 

5.2 Bilocal Moves: The Slithering-Snake and the Extended Reptation Algorithms 

A bilocal A^-conserving move consists in altering the configuration of two small groups of 
consecutive monomers. The groups are usually far from one another along the backbone of 
the chain. Typical examples are the slithering-snake and the extended reptation algorithm: 

• The slithering-snake (or reptation) algorithm removes a bond from one chain end. 
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adds a new one to the other end and shifts the inner monomers one bond up along 
the chain in direction of the new bond (Fig. 9a). As the positions of the inner 
monomers remain unchanged, the chains "sUthers" along its contour during the MC 

move (whence the name of the algorithm).^ 

• The extended reptation algorithm transports a kink or an end group via a slithering 
motion along the chain." Commonly utilized moves are: (1) "Kink-kink reptation", 
which deletes a kink at some position along the chain and inserts a new one at another 
position (Fig. 9b). (2) "Kink-end reptation", which removes a kink somewhere along 
the chain and adds two new bonds at one of the chain ends (Fig. 9c^). (3) "End-kink 
reptation", the inverse of "kink-end reptation" (Fig. 9c<— ). 

In the remainder of this section we will concentrate on the slithering- snake algorithm. 
Extended reptation is only discussed in comparison to the sUthering- snake algorithm. 

Implementation and Ergodicity. The slithering- snake and the extended reptation algo- 
rithms can be implemented in two ways: in a shrinkage-growth or a growth- shrinkage 
fashion. As growth- shrinkage is just the inverse of shrinkage-growth, we illustrate the 
procedure for the latter via the example of an isolated chain.'' For the slithering-snake 
algorithm one chain end is selected at random, the bond to its neighbor is cleaved, and a 
randomly chosen new bond vector is attached at the other end. If this move respects the 
excluded- volume condition in the athermal case and additionally passes the Metropolis test 
in the thermal case, it is accepted. Otherwise the old configuration is recounted. For the 
extended reptation algorithm the procedure is more complicated. Details may be found 
for the SAW on a hypercubic lattice in Refs. 51, 76 and for a continuum bond- fluctuation 
model in Ref. 77. 

Usually, shrinkage-growth is preferred to the growth- shrinkage procedure because it is 
computationally more efficient. The reason for this is illustrated in Fig. 10. The nested 
configuration of Fig. 10a would be frozen, if a new bond had to be appended before an 
end bond may be removed. However, it can be unraveled when shrinkage is attempted 
first. Thus, the shrinkage-growth algorithm is less plagued by -though not exempt of- 
non-ergodicity effects. An example is provided by the double cul-de-sac configuration of 
Fig. lOb.'^ It is frozen in the shrinkage-growth procedure, but not for the kink-end reptation 
move shown in Fig. lOc. In fact, kink-end/end-kink moves are known to be ergodic^^ (as 
well as other bilocal algorithms; see Ref. 50 for a thorough discussion). 

Should one thus abandon the sUthering- snake algorithm in favor of extended reptation? 
Usually, the answer is "No". For the SAW on the hypercubic lattice problems with ergod- 
icity arise to the constraints imposed by the small coordination number of the lattice. If 

^The slithering- snake algorithm was invented by Kron in the 1960's and later independently by Wall and Man- 
del.^'* For an overview of applications to SAW's see e.g. Ref. 43 and to off-lattice models see e.g. Ref. 75. 
"This generalization of the shthering-snake algorithm was first discussed in detail by Reiter.'* More recenfly, 
algorithmic and statistical properties of extended reptation moves were analyzed and their implementation was 
discussed in Refs. 50, 51. 

''For the multi-chain system the only difference to the isolated chain is that additionally one chain out of the n 
chains in the systems has to be chosen at random. 

°Non-ergodicity effects are less severe for the sUthering-snake algorithm than for the AT-conserving local moves 
discussed before. For the sUthering-snake algorithm the ergodicity class of a straight rod contains at least a 
fraction of Ar~(T~iV2 gf all SAW configurations, whereas this fraction is exponentially small for the local 
algorithms.^*'^' 
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Figure 10. Configurations of 2D SAW's to illustrate the ergodicity problem of the slithering-snake algorithm 
(a,b) and its solution via extended reptation moves (c). Panel (a) shows a configuration that cannot be moved 
by slithering-snake moves, if chain growth is attempted first. However, it is not blocked in the shrinkage-growth 
scheme. By contrast, the configuration of panel (b) is frozen for both growth-shrinkage and shrinkage-growth 
moves. Panel (c) shows that this configuration may be dissolved by extended reptation moves, e.g. by kink-end 
reptation if the chain end, where the kink is, happens to be selected for the attachment of the two bonds. 



many more bond vectors are a priori possible, as for the bond-fluctuation model or for 
(typical) continuum models, non-ergodicity should not represent a problem.'* 

Relaxation Time and Computational Complexity: Isolated Chains. One expects that the 
slithering-snake algorithm is able to decorrelate configurations more efficiently than a local 
updating scheme, the speed-up factor being roughly of order N. This expectation results 
from the following heuristic argument: The elementary move of the algorithm may be 
interpreted as a shift of all monomers along the contour of the chain. For the CM this 
curvilinear motion has two consequences: (1) The curvilinear diffusion coefficient 
should not depend on N, since all monomers are always shifted at once, irrespective of 
chain length. Thereby, the slithering-snake algorithm gains a factor of N in regard to the 
physical reptation dynamics, in which the curvilinear displacement is Rouse-like (Sec. 2.3). 
(2) An elementary move displaces the CM by ^ 6 along the chain backbone. After N such 
moves, the CM has diffused curvilinearly a distance of the order of the contour length 
L oc Nb. Thus, the relaxation time tn should be given by 

TAT = — =^ rjv - A^*^ (40) 



■^See pp. 283/284 of Ref. 43 for further discussion of that point. Contrary to SAW's, the equilibrium configura- 
tions of collapsed chains are typically (very) dense. Quite generally with increasing density, the slithering-snake 
or the extended reptation algorithm become less efficient, as the "free volume" to add new bonds decreases (see 
e.g. Ref. 78 for a comparison of various algorithms to simulate high-density polymer systems and the subsequent 
discussion). However, the recent study of Ref. 51 for 2D SAW's with N < 3200 at the 0-point suggests that 
extended reptation is almost as efficient as for pure SAW's with no attractive interactions. 
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Figure 11. Relaxation time tn versus N. tm is defined by gz{TN) = Rg- In the dilute limit (4> = 0.0078) 
Tjv ~ N'^^, as expected from Eq. (40). In the melt (</> = 0.5), the increase of Tj^r with N is stronger. The 
stretched exponentials are motivated by the activated reptation hypothesis:^^'*^'*' tjv ^ exp{0.8Ari/3) 
(bold line) provides a better description than tn «i exp(0.0747V-^/^) (dashed grey line). Adapted from 
Ref. 80. 



With respect to the computational complexity (page 22) one expects t^c « t^. There is 
no extra factor of N, as in the case of local moves, for the slithering- snake algorithm. The 
algorithm is bilocal. It takes a time of order 1 to check and update the chain ends.^ 

Note that Eq. (40) is independent of the conformational properties of the chain, contrary 
to Eq. (39) (which depends on i/). Thus, it should be valid for both 2D and 3D dilute 
polymer solutions as well as for dense melts. While for the slithering- snake algorithm*^ 
and for some extended reptation algorithms''-^*' the scaling found for tat is very close to 
Eq. (40), the behavior in dense systems is quite different.^^'*^" The influence of density on 
the sUthering- snake dynamics-^ has recentiy been studied by the bond-fluctuation model. 
The foUowing paragraphs briefly sunmiarize some results of this work. 

From Dilute Solutions to Dense Melts: A Case Study by the BFM. Reference 80 describes 
simulations for athermal systems containing chains of length 16 < < 1024 at different 
volume fractions (j). (j) ranges from dilute solutions to dense melts {(j) « 0.5; see footnote 
on page 22). Figure 11 compares the relaxation time rjv in dilute solution with that in 
the melt. In dilute solution, the simulation results agree with the prediction of Eq. (40), 
Tjv iV~^. This implies that the assumption of independent, free diffusive motion, which 
underlies Eq. (40), is well borne out. If this assumption was also true in the melt, the 
sole effect of density would be to slow down the monomer mobility m. However, the 



'^Here, we assume that the time to shift the monomer index along the chain is implemented in a way that it also 
takes a time of order 1 only. 

■'^For local moves the set of allowed bond vectors automatically prevents bonds from crossing each other in the 
course of the simulation (see Sec. 4.1). If sUthering-snake moves are considered, the uncrossabiUty of the bonds 
has to be checked explicitly to avoid configurations which caimot be attained or unraveled by local updates. Bond 
crossing can occur if vectors from the classes [2, 2, 1] or [3, 1, 0] are selected.^' 
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dependence of rjv on should not change. Figure 1 1 shows that this is not true. At 
= 0.5, Tjsf increases exponentially with N. This strong slowing-down of the dynamics 
reflects correlations between the motion of the chains. 

The importance of such intermolecular interactions for the polymer dynamics was first 
discussed by Deutsch.'^ However, Deutsch goes beyond a mere interpretation of the dy- 
namic properties of the slithering-snake algorithm. He identifies the slithering dynamics 
with the physical dynamics along the primitive path in the reptation model (see Fig. 5). 
This suggests an attractive application: The sUthering- snake algorithm mimics the back 
and forth reptation motion of real chains without modeling the (time consuming) local 
monomer fluctuations around the primitive path. It focuses on the long time behavior of 
very large chains, where all of these local motions have already relaxed. This suggests 
that the slithering-snake dynamics may be interpreted in terms of theories proposed for the 
dynamics of strongly entangled polymer melts, such as the one of Deutsch.^' 

The main results of this theory may be summarized as follows: A chain can reptate 
through the network of its neighbors only as long as the end monomer does not enter a 
dense region which prohibits any further forward move. The only way out of the trap is to 
partially retract and to explore the environment for new pathways. These intermolecular 
interactions create a free energy barrier which temporarily localizes the chain in the region 
it initially occupied, and protracts the relaxation. Further relaxation in a dense region 
could only occur if the chain end encounters another end which moves out of its way. 
This implies that the portion of the chain, which altered its initial configuration while 
exploring the environment, should span the typical distance between chain ends ciend- Let 
there be g monomers in this portion. Then, by exploiting the ideality of the chains in 
the melt, we have g = (dend/^)^ ^ A^2/3 because the density of chain ends scales as 
p/N (X Thus, g is large for long chains. If we now assume that the monomers 

have to overcome the free energy barrier gA/j,, where A/z is the difference in the monomer 
chemical potential between the newly explored environment and the region of the initial 
chain configuration, the barrier is large and the relaxation dynamics should be activated. 
Thus, tn oc N'^ exp[const A/"^/^]. This is the main prediction by Deutsch. The assumption 
of a finite A/7, was challenged by Semenov^^ who suggested that the barrier is due to 
fluctuations of the molecular field rather than to a permanent chemical potential difference 
(see also Ref. 83). This picture impUes that the barrier should be proportional to ^ so 
that TAT cx iV^ cxp[constA^i/^]. 

The simulation data of Fig. 1 1 appear to agree with the latter prediction better than with 
the original one of Ref. 79 (at least for the chain lengths simulated up to now). Certainly, 
more work is needed to test these predictions. 

Slithering-Snake versus Local Moves. From a merely computational point of view Fig. 1 1 
appears to indicate that the slithering- snake dynamics is not very efficient in equiUbrating 
dense melts. Its relaxation time increases with N more strongly than a power law which 
is typically found for local updating schemes.** Nevertheless, simulations of the BFM for 
short chains (N = 10, 20) suggest that the slithering- snake algorithm decorrelates melt 
configurations ((/) w 0.5) very efficiently.**'*'*'^ 

This point certainly needs more studies. Work in this direction was done in Ref. 80. 
Figure 12 shows a prelinninary result for the diffusion coefficient I? jv as a function of chain 
length. Dn was obtained from simulations employing a mixture of local and slithering- 
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Figure 12. Spatial diffusion coefficient Djv versus A'^ for the BFM at = 0.5. Different ratios ui of local to 
slithering-snake moves are compared: uj = corresponds to pure slithering-snake dynamics, ui = ooto the pure 
local dynamics. The diffusion coefficient is scaled by N/A, where A is the acceptance rate. For both local and 
slithering-snake dynamics the acceptance rate is roughly A x 0.1 for all Na.t<p = 0.5. The data for a; <^ 1 and 
for 2> 1 are very similar to the pure slithering- snake {lu = 0) and the pure local limit (ui = oo), respectively. 
For f5i 8, A^Djv /A is approximately independent of A^. This may define a reasonable choice of lu for efficient 
equilibration of longer chains by local and slitheiing-snake moves. Adapted from Ref. 80. 

snake moves. This introduces, as an additional parameter, the ratio uj of local to slithering- 
snake moves. The figure indicates that pure slithering-snake dynamics (uj = 0) equilibrates 
short chains more efficiently than pure local dynamics (uj = oo), in accord with the obser- 
vations made in Refs. 84,85. By contrast, with increasing chain length slows down ex- 
ponentially for the slithering-snake algorithm, as expected due to 13 at ~ R^/tn, whereas 
the local dynamics exhibits a crossover from Rouse-like, ^ 1/iV, to reptation-like 
behavior. Dm ^ 1/N~^-* (see Sec. 2.3). If this trend persists, the pure slithering-snake 
algorithm will become inefficient to equilibrate long chains. However, one can speculate 
that the addition of local moves weakens the confinement imposed by neighboring chains 
on the slithering-snake dynamics in the large-A^ limit. Indeed, this seems to be borne out 
by the data. For short chains (N < 64) Dn decreases monotonously with increasing oj, 
since local moves are less efficient in exploring the configuration space and the confine- 
ment is negligible. As N increases, one finds a non-monotonous behavior. The dynamics 
first becomes more rapid, as local moves are added. This effect appears to saturate at 
u! w 10. Larger values of lu causes the diffusion coefficient to decrease again strongly (at 
fixed N). This implies that a judicious (model-dependent) choice of uj is crucial if one 
wants to equilibrate a melt of long chains efficiently by mixing local and slithering-snake 
moves. 

Remark. The efficiency of the slithering-snake algorithm (with or without local moves) 
or of the extended reptation algorithm deteriorates considerably as (f> approaches 1, since 
there is not sufficient space for the growth step. If one is interested in these high den- 
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sities,^ an alternative simulation method may be provided by reptation moves including 
a "walker". ^^'^^ A "walker" is defined as an isolated monomer (or as a small group of 
monomers). In the MC move, the "walker" attaches to a chain in its neighborhood, which 
then releases a monomer somewhere along its backbone, yielding a new "walker" at a 
different position than the original one. Since the walker can be created by cleaving a 
monomer from a chain, the algorithm works even at (f> = 1. In order to preserve monodis- 
persity the update of the configuration is finished if the "walker" attaches again to the chain 
it was originally cleaved from. 

5.3 Non-Local Moves: The Pivot Algorithm 

A non-local A^-conserving move attempts to update a chain portion of order N. If success- 
ful, it drastically modifies the chain dimension. Global properties, such as the end-to-end 
distance, should therefore relax within a few (A^-independent) steps so that the critical 
slowing-down is largely moderated -if not removed.'* 

This appealing feature makes the search for appropriate non-local moves very attrac- 
tive. However, not every conceivable move turns out to be efficient. This is mainly due to 
two reasons: 

1. A drastic modification of the chain configuration is much more likely to violate the 

excluded-volume constraint than local or bilocal moves do. So, we expect the ac- 
ceptance rate of non-local updates to decrease with N. The challenge consists in 
inventing moves whose acceptance rate does not rapidly vanish as A^ — > oo. 

2. A non-local move typically requires a CPU time of order N (checking self-avoidance, 
updating the configuration if accepted), in contrast to order 1 for a local or a bilocal 
move. The extra factor N must be compensated by a very efficient decorrelation of 
the configurations to justify the use of non-local moves. 

A paradigm for a non-local algorithm, which satisfies these criteria, is the "pivot algo- 
rithm".' 

Pivot Algorithm for the SAW. The elementary move of the pivot algorithm works as fol- 
lows (Fig. 13): First, one randomly chooses a monomer i and a symmetry operation (e.g. 
a rotation, a reflection, etc.). The monomer serves as a "pivot point" for the symmetry 
operation which turns the chain portion comprising the monomers i + 1, . . . , A to a new 
position, while the other piece of the chain (monomers 1, . . . , i) remains unchanged. The 

^Note that </> 1 means that the polymer melt has zero compressibility kj-. By contrast, real polymer melts are 
compressible, with k^TpKT being of the order 10^^ for temperatures above the glass transition or crystallization 
temperatures.'*'' For the BFM^^ at 0.5, kgTpKx ~ 0.2 and for the bead-spring models of Sec. 4.2 at 
p Si 0.9, 10^^ < k^TpKT < 10~l (see Refs. 24, 87). Thus, it appears that the limit — > 1 is not needed to 
model dense polymer melts. 

''This idea is related to that of cluster algorithms'^' employed in MC simulations of spin systems near criticality. 
Close to the phase transition, the spins are strongly correlated. They form cluster of size ^ (= correlation length, 
corresponding to i?e in the polymer problem, see Sec. 1). A cluster algorithm finds one cluster (Wolff algorithm) 
or all of them (Swendsen-Wang algorithm) and updates all spins of the cluster(s) at once. This strongly reduces 
or even eliminates, in favorable cases, the critical slowing-down.'^'*^ 

*The pivot algorithm was invented by Lai in 1969. A comprehensive discussion of the algorithm may be found 
in Refs. 34, 89. 
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Figure 13. Illustration of the pivot algorithm. A monomer i (= pivot point) is chosen at random. It divides the 
chain into two pieces: The monomers 1, . . . ,i remain fixed, while the monomers i + . . . , N are translated to a 
new position via a randomly chosen symmetry operation (rotation, reflection, etc.). In the example of the figure, 
a 180° rotation around the monomer i is shown. 



proposed move is accepted if the resulting configuration is self-avoiding. Otherwise, it is 
rejected and the old configuration is recounted. 

Qualitatively, the pivot move resembles an attempt to construct a SAW of ^^A^ mono- 
mers by joining two SAW's of monomers at the pivot point. The probability for 
the resuh to be self-avoiding should scale as Z{N)/ Z'^{N/2) - Af^t'^^i) [see Eq. (8)]. 
This heuristic argument suggests that the acceptance rate vanishes as 7V~~" in 2D 
= 43/32)^' and as iV«-o.i58 jjj ^ 1.158)^'. Although these estimates are quan- 
titatively not very accurate, they correctly predict the qualitative trends: The acceptance 
rate decreases with increasing chain length as a power law A^^" and the exponent a is 
larger in 2D than in 3D (2D: a « 0.19, 3D: a « 0.11).^'' 

Relaxation Time and Computational Complexity. Fortunately, the numerical value of a is 
small, implying that even for long chains, e.g. A^ — 10^, every A^^-^^ « 3.5th move is 
accepted. Since a successful move implies a huge modification of the conformation, one 
can expect global properties to relax after a few steps. So, the relaxation time ta? scales as 
tn ^ N". This increase is distinctly slower than that of all algorithms discussed so far. 
Due to this property and due to the fact that the pivot algorithm is known to be ergodic**^ it 
has become very popular (see e.g. the compilation of references in Ref. 90). Currently, the 
pivot algorithm is considered to be the most efficient algorithm for studying configurational 
properties^ of isolated SAW's.^i'**'' ''^ 

The pivot algorithm quickly decorrelates global quantities, such as the end-to-end dis- 
tance. However, it is not as efficient for local properties: The conformation of a specific 
monomer is only altered if this monomer is selected as a pivot point and if the move is 
successful. A successful moves takes a time of order A^" and, as the chain consists of A^ 
monomers, the decorrelation time of a local observable should scale as rioc ~ A^^+". This 

^By configurational properties we mean quantities characterizing the chain dimension, such as /?c, Rg, etc. A 
very accurate estimate of the critical exponent u in 3D (u = 0.5877 it 0.0006) was obtained by the pivot 
algorithm." By contrast, it appears difficult to measure precisely the partition function Z(N), and so the ex- 
ponent 7 with the pivot algorithm. For this purpose, other algorithms, such as the "join-and-cut" algorithm^' or 
chain-growth algorithms'^, are better suited. The current best estimate for 7 in 3D is 7 = 0.1575 it 0.0006."' 
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extra factor N is felt if one starts from an arbitrary initial configuration. Full equilibration 
on all length scales is required before large-scale equilibrium properties may be sampled. 
The equilibration time must be longer than the longest relaxation time in the system, i.e, 
than rioc. 

For a non-local algorithm the computational complexity is a particularly important 

quantity because inefficient implementations may ruin the advantage gained by fast decor- 
relation. The most naive check for self-avoidance would take a time of order N"^ so that 
Tec = N^tn ~ N'^~^", comparable to the slithering-snake algorithm [Eq. (40)]. Obvi- 
ously, a faster check is called for. In Ref. 89 it was argued that, by starting at the pivot 
point and working outwards, self-intersections may be detected in a time of order iV^~". 
This procedure must be repeated ~ N" times to obtain one accepted pivot. So, the time 
required per accepted pivot scales as ^ N.'^ Once the pivot is accepted, we still have to 
update the monomer positions which also takes a time of order N. So, in total we find 
Tec ~ Ntn ~ N^'^" for global properties and rl"^ ~ A^'rioc ^ N'^~^°' for local proper- 
ties. The estimate for t'°'^ is again comparable to the slithering-snake algorithm [Eq. (40)]. 
Therefore, one could also use shthering- snake moves to engender an initial, equihbrated 
configuration for the pivot algorithm. 

5.4 Non-Local Moves in the Melt: The Double-Pivot Algorithm 

Due to its efficiency in decorrelating SAW configurations it is tempting to apply the pivot 
algorithm also to other situations, such as the collapse transition, SAW's in confined geom- 
etry, or dense polymer melts. However, in these cases, the algorithm becomes inefficient. 
The non-local moves either lead to large energy differences (collapse transition) or violate 
the excluded-volume condition (SAW's in confined geometry, dense melts) so that they are 
rejected. 

Should one therefore give up the idea of using pivot-like moves, say, for dense melts? 
Recent work suggests that this conclusion might be wrong. Instead of pivoting a piece of 
one chain to a new position the MC move can involve two chains. Such a move was termed 
"double-pivot (DP)" algorithm.' The basic idea of the algorithm is to cleave simultaneously 
a bond in a chain and in one of its neighbor chains, and to reconnect the monomers such 
that the chains remain monodisperse.™ The algorithm works as follows (Fig. 14): 

1 . A monomer, say monomer i in chain a, is chosen at random. Around this monomer the 
neighborhood is inspected to find bridgeable neighbors on other chains b. A bridge- 
able neighbor is defined by the following requirements: 

(a) It must be possible to connect the neighbor to i by a bond vector. This imposes a 
restriction on the intermolecular distance between i and its neighbor. In the ex- 

'=The pivot algorithm may be implemented so that the time reqviired to obtain an accepted move is of order AT^ 
with (J < 1.* 

'The double-pivot algorithm has been proposed recently in Ref. 93. In part, this vifork was motivated by a novel 
chain-bridging algorithm which was successfully employed in atomistic simulations of long polyethylene chains 
(see Ref. 94). Our presentation of the DP algorithm is inspired by the discussion of Ref. 94. Due to the newness 
of the algorithm the implementation that we propose might turn out not to be the most efficient one. 
"In general, connectivity-altering moves between arbitrary monomer pairs of neighboring chains lead to a dis- 
tribution of chain lengths, i.e., to polydispersity. The width of the distribution can be controlled by introducing a 
chemical potential. See Sec. 3.3 of Ref. 43 for further discussion in the context of lattice models and e.g. Ref. 94 
for an appUcation to atomistic MC simulations of polyethylene. 
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Figure 14. Illustration of a double-pivot (DP) move on a square lattice. A DP move flips bonds between two adja- 
cent chains a and b. To this end, the neighborhood of monomer i in chain a is inspected to find potential bridging 
sites on chain b. To preserve monodispersity a potential bridging site has to satisfy the following conditions: (a) 
The {i + l)th monomer of chain b, denoted by i' + 1 in the figure, must be separated by a distance of the bond 
length from monomer i. (b) The same condition must also hold for the distance of the ith monomer of chain b, 
denoted by i' , from monomer i -f 1 of chain a. Between these four monomers a connectivity-altering move is 
attempted. The bonds from t to j -f 1 and from i' to i' + 1 are broken, and new bonds, i with t' + 1 and i' with 
i + 1, are created. The proposed move is accepted or rejected according to Eq. (43). From the figure it is clear 
that the DP algorithm can only be carried out if there are matching monomers on a neighbor chain within the 
distance of a bond. This is the more likely, the higher the concentration of the solution. Therefore, the algorithm 
(presumably) works best in concentrated solutions or melts. If successful, the move entails a drastic change of 
the chain configuration. 



ample of Fig. 14 it must coincide with the lattice constant. For the BFM it must 
be among the set of allowed bond vectors [Eq. (33)], whereas, for a continuum 
model, the bond energy resulting from taking the intermolecular distance as a 
bond vector should not be so large that the proposed bond would never occur in 
equilibrium. In the latter case, it might be necessary to reduce the strength of 
the bond potential, e.g., the force constant of Eq. (34).^-' 

(b) To maintain monodispersity the neighbor must be either the (i — l)th or the 
(i + l)th monomer of the chain h. We distinguish the monomers of chain 6 by a 
"prime", e.g. i' , from those of chain a. 

(c) If it is monomer i' ± 1, monomer i' must be separated from monomer i ± 1 of 
chain a by a distance which satisfies condition (a). 

Using these three criteria we determine the total number of bridgeable neighbors of 
monomer i, Nup{i,x'). If Nr,p{i,x') — for all i, the configuration x' must be 
updated by local (or bilocal) moves to bring the monomers in more favorable positions 
for bridging. 

2. A double-pivot move is initiated by randomly selecting one of the bridgeable neigh- 
bors, say i' + 1 in chain b, from A^dp(*, x'). Then, the bonds between i and i + 1 and 
between i' and i' + 1 are cleaved, and one attempts to create new bonds between i 
and i' + 1 and between i' and i + 1. This move just switches the connectivity while 
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preserving the chain length, and is proposed with probability 

1 



Pj^ix' ^x) = ^r-7T^ ■ (41) 



P,ro{x ^ x') = , , ^„ (42) 



3. To satisfy detailed balance we have to determine Ppmix x') of the reverse move. 
As the forward move only alters the connectivity between the chains, but does not 
displace the monomers, the number of bridgeable neighbors of a specific monomer 
remains unchanged. That is, A^dp(«, x) = -/Vdp(z, x'). To reverse the forward move, 
we have to select monomer i' + 1 on chain b and its bridgeable neighbor i on chain a. 
This occurs with probability 

1 

iVDp(i'"+l,a;') 
so that the Metropolis criterion reads 

acc,.'-«, = ..(,^,;^H<^e-l<'<-.-<-'.l). (43, 

The difference U{x) — U{x') is the local change in energy due to the switching of 
the bonds between chains a and b. 

The steps 1.-3. may be repeated several times. However, the number of iterations should 
not be too large. Otherwise it is likely that an accepted move annihilates one of its pre- 
decessors by performing the transition between two chains in the reverse direction. To 
avoid this inefficiency it is important to mix up the local configuration of the system. This 
may be achieved by e.g. local MC moves or by combining'^^ the DP algorithm with MD 
simulations. 



6 Monte Carlo Methods for Polymers: Rosenbluth Sampling and Its 
Modem Variants 

The first MC method to simulate a SAW was "simple sampling"."'^ This static method 
(Sec. 3) works as follows: 

1 . Place the first monomer at the origin, randomly choose a bond vector, and append it 
to the monomer. 

2. Choose the next bond vector, again randomly, connect it to the second monomer, and 
check the self-avoidance (Fig. 15). 

3. If the chain is self-avoiding, the random growth process may be continued. If not, the 
self-avoiding piece of the SAW, obtained up to this point, must be discarded, and we 
have to start from scratch at the first step again. 

4. The steps 1. to 3. are repeated until a SAW of the desired length N is obtained. Then, 
data analysis may be done. 

5. Repeat steps 1. to 4. to gather sufficient statistics. 
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Figure 15. Illustration of the simple-sampling (a) and the Rosenbluth-Rosenbluth (RR) methods [(b) and (c)] on 
a square lattice. The coordination number of the lattice z (= 4) defines the number of possible bond vectors 
{b}. In simple sampling, all of these vectors have the same a priori probability. Thus, the bond vector 
(= — ri_i) from monomer i — 1 to the new monomer i can also point in the direction of already occupied 
lattice sites. This leads to the attrition problem. The RR method strongly reduces the attrition by taking 
only from the open directions ki (= 2 in the example of the figure). The new bond vector can be chosen either 
uniformly from the ki possibilities, as indicated by the dashed lines and dashed circles in panel (b), or according 
to the local Boltzmann factor [Eq. (48)], if e.g. an attractive monomer-monomer interaction is present [filled grey 
circle in panel (c)]. 



Apparently, completed SAW's are independent of one another and occur with the same 
probability Ps{x) = Psavj{N). This is the main advantage of a static MC method. The 
main disadvantage of simple sampling is that Psaw(^) becomes exponentially small for 
large N. To see that, let us calculate Psaw(^) for a SAW on a hypercubic lattice. The 
hypercubic lattice has the coordination number z = 2d. Thus, the number of random 
walks (RW's) starting at the origin and having iV — 1 steps, is Zrw — z^^^. Zrw is 
the partition function of the RW. Out of these z^^^ random-walk configurations simple 
sampling selects those which are self-avoiding. As there are Z (i^ N^~^, Eq. (8)] such 
configurations, PsnviiN) is given by 

Psaw(A^) ^ ( ^ ) = N^'^ e-^^ - e-^^ {N large) , (44) 

where A = ln(z//i) is called "attrition constant". The attrition constant is A > (in 2D 
and 3D)," reflecting that the monomer partition function, i.e., the number of ways to place 
a monomer on the lattice, of a RW (= z) is, due to the neglect of self-avoidance, larger 
than that of a SAW (= fj,). 

Equation (44) illustrates that simple sampling is not an efficient simulation method 
for SAW's. On average, e'^^ random walks have to be constructed to obtain one SAW 
("attrition problem"). As A = 0.416 and 0.248 for the 2D and 3D hypercubic lattices, this 
number of constructions is prohibitively large already for N > 50. Even if one modifies 
the construction method by avoiding immediate backfolds, thereby replacing z by z — 1 
("Non-Reversal Random Walk (NRRW)"), the exponential attrition remains. Generation 
of SAW's with N > 10^ is still unfeasible. Therefore, all alternative simulation techniques 
have to alleviate this attrition problem. 

"For a compilation of attrition constants see Refs. 34,43. If d — > oo, A goes to zero as A — > l/2d for the 
hypercubic lattice. 
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6.1 Inversely Restricted Sampling: Rosenbluth-Rosenbluth Method 

The attrition problem arises because simple sampling chooses blindly from the nearest- 
neighbor sites to place a new monomer. A more clever algorithm could scan the local 
environment around the last monomer and exclude those trial directions which lead to self- 
intersections. The position of the next monomer i can then be chosen with equal proba- 
biUty from the remaining fc, open directions (Fig. 15).° This method, known as "inversely 
restricted sampling" or "Rosenbluth-Rosenbluth (RR) algorithm",^' strongly reduces the 
attrition^ at the expense of introducing a bias. A SAW is not generated with uniform prob- 
ability, but with probabiUty 

^ 1 

^^('^) = nTr (^1 = 1)- (45) 

Equation (45) shows that configurations with small fc^'s have a higher probability of occur- 
ring. This bias toward dense configurations in the production of a SAW must be corrected 
in its analysis by the weight W{x) oc 1/P^{x) when calculating observables [see Eq. (18)]. 

The RR algorithm is a static MC method. As such, it has the advantage that suc- 
cessively generated SAW's are independent of each other. All problems of decorrelating 
configurations, discussed in Sec. 5, are absent by construction. On the other hand, Eq. (45) 
also points to the major difficulty of the method. The RR method favors dense config- 
urations which are not representative of long SAW's. Thus, Ps{x) differs from Peq(a^)- 
As N increases, the difference becomes more pronounced. To compensate the discrep- 
ancy between the two probabilities the distribution of weights must become broad: Dense 
configurations have small weights and open chains have in general larger weights. A de- 
tailed analysis of this problem was undertaken by Batoulis and Kremer^^ They showed 
that the distribution of weights, obtained from AI repetitions of the RR method ({ ^^(Xtti), 
m = 1, . . . , M), is dominated by few configurations having the largest weights. The most 
relevant SAW configurations have, however, smaller weights. To sample this portion of 
the weight distribution sufficiently, M has to become very large [see the discussion of 
Eq. (22)]. This problem makes the RR algorithm not suitable for the simulation of long 
SAW's."" 

However, the RR algorithm should be well suited if the bias introduced by the sampling 
engenders configurations which are close to the physical ones. That is, if the equilibrium 
configurations are less swollen than those of a SAW. As this is the case close to ©-point in 
3D (see Sec. 1), one might expect the RR algorithm to be more efficient for T ps Tq. In 
fact, this expectation is nicely borne out. The biased sampling of the RR method produces 

°Here, we assume < fc^ < 2: — 1. If all neighbor sites are blocked (ki = 0), the configuration obtained until 
then is trapped. It must be discarded, and the construction resumes from the beginning. The first monomer can be 
placed anywhere on the lattice. If the construction always starts, say, from the origin, we set fci = 1 in Eq. (45). 
''See Ref. 95 for a detailed statistical analysis of the RR algorithm'* in the context of simulating SAW's. Chap- 
ter 1 1 of Ref. 3 1 gives a discussion in the larger context of free energy calculations with appUcations to discrete 
and continuous chain models. 

^In the RR algorithm, attrition occurs because the growing chain may be trapped, i.e., fe, = for some i. In 
practice, this does not appear to be a serious problem in 3D for high-coordination lattices.'^''^ It was suggested 
that effects of trapping should become visible only for N > 10'*.'^ This implies that there is stiU exponential 
attrition, as AT -+ 00, but at a much lower rate than in Eq. (44). 

''In the simulations on the (high-coordination) FCC lattice of Ref. 95 the systematic error due to the weights 
rendered a precise determination of Rg impossible for N > 200. 
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an effective atoaction between the monomers which closely resembles that of 6-chains. 
An important consequence of this special property is that the weights are nearly compen- 
sated by the Boltzmann factor. Thus, the RR method was employed to study properties 
near the 0-point (see e.g. Refs. 10,11 of Ref. 97), and it also represents the core of a 
modem algorithm, the "Pruned-Enriched Rosenbluth Method (PERM)''''". 



6.2 Pruned-Enriched Rosenbluth Method (PERM) 



The simulation of a polymer chain close to the 9-point requires the introduction of an 
attraction between the monomers to compensate their mutual repulsion. Typically, these 
thermal interactions are modeled by a short range inter-monomer potential.* As a chain is 
grown according to the Rosenbluth scheme, the presence of the potential implies that the 
internal energy of a chain changes: 

Ui{bi-i) = U{ri, . . . ,fi-i +bi-i) - U{fi, . . . ,fi-i) [ui{bo) := ?7(ri)] . (46) 

Here, U{fi, . . . , ri-i) is the potential energy of chain having (i— 1) monomers and bi-i 
is the bond vector from the (i — l)th to the ith monomer. A priori, there are different ways 
to incorporate Ui{bi^i) in the construction. Two possible choices are the following:"*^ 

1 . The first method is the classical Rosenbluth scheme. The position of the ith monomer 
is chosen from the free neighbors with uniform probability [see Eq. (45)]. Thus, the 
weight of the new chain configuration x {= fi, . . . , rjv) is given by * 



W{x) 



N 

n 



^ g-/3ui(6i_i) 



(47) 



2. An alternative consists in including the Boltzmann factor in the probabiUty P^^i for 
placing the iXh monomer. Let {6} denote the ensemble of possible bond vectors. For 

the hypercubic lattice, {6} coincides with the number of lattice directions z, for the 
BFM it is given by Eq. (33). Then, we may write for P^ i (Fig. 15) 
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where the normahzation Wj reduces to fcj if only excluded-volume interactions are 
taken into account (SAW Umit). This implies that the weight W{x) is given by 



W{x) = 



.-0U{x) 

'pM' 



N 
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-/3«,(&i_i)- 



N 

i=l 



(49) 



"For the SAW on the hypercubic lattice the attraction is usually implemented between non-bonded nearest neigh- 
bors.''^ In a simulation of the ©-transition with the BFM, a square-well potential of range \/6 (in units of the lat- 
tice constant) was used.** This choice ensures that the first peak of the pair-distribution function is encompassed 
by the range of the potential. A more compUcated choice was made in Ref. 98. In the continuiun, a Gaussian 
chain model with a non-truncated LJ-interaction was extensively studied (see Refs. 64,99 and references therein). 
*Note that, contrary to the definition of W{x) in Sec. 3, Eq. (47) does not include the factor of the (vmknown) 
partition function. 



36 



Both methods have been used to simulate 0-polymers via the Rosenbluth algorithm. ^ 

Despite the fact that, precisely at Te, the RR configurations more or less coincide 
with the equilibrium configurations, the accuracy of the method deteriorates for TV > TVc." 
Again, the reason is that Ps{x) does not perfectly agree with Pgq (a; ) . As a result, the weight 
distribution becomes so broad that chains with the biggest weights dominate the sample 
(see Fig. 2 of Ref. 47). Due to Eq. (18) this leads to a large variance of the computed 
observables. 

To improve the accuracy one has to reduce the variance. Grassberger proposed a clever 
way to achieve this."*^ Assume that we have constructed a chain up to monomer i (1 < 
i < N) via the RR method. This chain has the weight Wj which we want to prevent from 
fluctuating too much. That is, if Wi exceeds a lower or an upper bound, we interfere in the 
following way: 

1. If Wj < Wf, we "prune" the sample: If Wi decreases below the threshold W^", a 
random number < C < 1 is uniformly drawn. If C < 1/2, the chain is removed. 
Otherwise, it is kept, its weight is doubled (Wj 2Wi), and the step-by-step growth 
continues. 

2. If Wi > W^, we "enrich"" the sample: If Wi exceeds the upper bound W^^, c copies, 
typically 2,'*^" of the configuration are made, each of which is given the new weight 
Wi ^ Wi/ c. These copies are then grown independently of each other. 

This control of the weight distribution within the RR algorithm was termed "Pruned-En- 
riched Rosenbluth Method (PERM)".'*'^ 

Of course, the question arises of how to choose the bounds W^ . Here, it is important to 
note that neither the pruning nor the enrichment step introduces any bias. In the calculation 
of the sums in Eq. (18) the increase of Wi by pruning is compensated by the probability 
1/2 with which the configuration is retained, and the decrease of Wj in the enrichment is 
compensated by the number of copies c. Thus, we are free to choose the bounds W^. Bad 
choices can "only" render the method inefficient, but not incorrect. In order to determine 
optimum values for W^ the following procedure was proposed (for temperatures that are 
nottoolow):^'''!'"' 

• First, one chooses H'^ = and W^ very large. That is, one performs a simulation 
via the original RR method. This simulation yields the weights Wi for i = 1, . . . , N. 
First estimates for the bounds W^ are then determined by Wf = Wi and W^ = 
C+Wi with C+/C- « O(l)-O(10). 

• These estimates are refined "on the fiy". Imagine that we have obtained Mi config- 
urations of chain length i from the simulation. Then, we first calculate the partition 

"The critical chain length Nc depends on the model. For a simple cubic lattice it is « 10^ .'^^ 
"Enrichment is a classical technique for simulating SAW's. Briefly, it works as follows: If a chain survives the 
s step, c copies of its configuration are made, which serve as independent starting points for further growth. The 
method may be implemented in a "breadth-first" or a "depth-first" fashion. The former impUes that all copies are 
first grown to size 2s before the entire sample is copied again. By contrast, the latter method tries to complete 
the construction of one copy up to chain length N before passing to next one. The pros and cons of the two 
implementations are discussed in the context of the PERM in Ref. 47. More details about enrichment may be 
found in Refs. 34,43. 
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function by [a; = (ri, . . . , fi), see Eq. (18)] 

^. = /d.e-*. . ± i: . ± ^ , (50, 

771—1 ■ ' 77X — 1 

and from that, we determine the new bounds by W^^^ = C^Z^. 

Applications of the PERM. The PERM was invented to simulate the transition from an 
excluded- volume to a collapsed chain at the 6-temperature Tq. Theoretically, the transi- 
tion is usually identified with a tricritical point in the limit N — > oo.^^ A tricritical point 
exhibits mean-field behavior in 3D. Thus, one expects that u = 1/2 and 7 = 1 [Eq. (8)] 
at Tq for N 00. This asymptotic large- A'^ behavior is supplemented by (universal) cor- 
rections of order 1/lnN for finite chain length (see e.g. Chap. 21 of Ref. 16 for a good 
discussion). A significant test of the theory therefore requires to study very long chains. 

In Ref. 47 such a test was attempted. Grassberger performed a comparative study of 
various models employed in the Uterature: a SAW on a simple cubic lattice with attrac- 
tive nearest-neighbor interactions, the BFM with two versions for the attractive monomer- 
monomer interactions,^"'^*^ and a LJ bead-spring model^"^. Using the PERM these models 
could be simulated with high precision and, partiy, with much longer chain lengths than 
studied before. These simulations yielded refined estimates of Tq for the various models, 
confirmed the mean-field-like asymptotic character of the 0-point, but also showed that 
the leading-order logarithmic corrections cannot explain the finite- A'' behavior found, even 
for N = 10^'" 

On the technical side, it was found that the selection of the position of the next 
monomer i with uniform probability from the ki open directions is only the best choice 
for the SAW on the simple cubic lattice (first method of page 36). As alluded to at the 
end of Sec. 6.1, this is due to a near cancellation of the Rosenbluth weight and the Boltz- 
mann factor: Many nearest-neigbor contacts lead to a low Rosenbluth weight, but to a large 
Boltzmann factor, and vice versa. For more long-range or more complicated interactions 
the degree of cancellation need not necessarily be the same. In fact, for the BFM it was 
found in Ref. 47 that a selection of the next monomer position according to the Boltzmann 
factor [Eq. (48)] is more efficient. Similar approaches were also used to study e.g. simple 
models of proteins. 

The preceding discussion appears to suggest that the PERM is a single-chain technique. 
This is not true. We just quote two recent examples. The PERM was utilized to simulate 
the denaturation transition of a simple model for double- stranded DNA (two SAW's).'*'^ A 
truly multi-chain system was studied in Ref. 103. This work is concerned with the phase 
diagram of semidilute polymer solutions for T < Q (see Fig. 3). For a review of these and 
other applications see Ref. 100. 

"These findings elicited furtlier tlieoredcal"" and numerical work."'"" On the theoretical side, subleading 
corrections of order lii(lii Af)/(ln Af)^ were calculated and found to be as large as the leading 1/lnN term. 
On the numerical side, MC simulations"" for N < 10* of a NRRW, including weakly attractive two-body, 
but repulsive three-body interactions, were performed. This model shows logarithmic corrections which are 
much weaker than those foimd in Ref. 47 and are roughly compatible with the theoretical predictions. However, 
Ref. 101 stresses a problem in the analysis of the ©-point. To estimate Tq, an infinite-AT property, precisely 
from the simulations one has to rely on the theoretical predictions for the finite- AT corrections to extrapolate to 
AT —» 00. To this end, the simulated chains must be long enough for these corrections to apply. This regime 
appears to be very hard to attain, even for N ~ 10®. 
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6.3 Conflgurational-Bias Monte Carlo and Recoil-Growth Algorithm 

An alternative multi-chain MC scheme, incorporating the RR method, is the Configura- 
tional-Bias Monte Carlo (CBMC) algorithm.^ Contrary to the PERM, CBMC builds up 
a new chain step-by-step without controlling the weights Wi "on the go". It is only after 
a successful construction that the resultant bias is removed: The new chain is accepted 
according to the Metropolis criterion with a probability dependent on the total Rosenbluth 
weights of the new and the old chain configurations. This additional test warrants sampling 
from the Boltzmann distribution. 

A recent extension of CBMC is the "Recoil-Growth (RG) algorithm".^ Contrary to the 
RR method, which only looks ahead one step while constructing a chain, the RG algorithm 
uses a more sophisticated growth procedure. It places a long retractable feeler at the head 
of a growing chain. The feeler spys out the environment to find favorable pathways for the 
chain construction. The efficiency of the method resides in the fact that the growth does not 
terminate if the feeler encounters a trap. It merely recoils back from the trap and pursues its 
search in a different direction. After the construction is completed, the new chain replaces 
the old one, just as in CBMC, with a probability determined by their respective weights 
according to the Metropolis criterion. 



6.3.1 Coniigurational-Bias Monte Carlo (CBMC) 

To illustrate the CBMC scheme in more detail we consider a solution of SAW's on a lattice 

(Fig. 16): 

1. Given the initial configuration x' of the system we randomly select a chain and one 
of its monomers. Let this be the monomer n' of chain 6 (1 < n' < A''). The confi- 
guration of the chain portion i = n' , . . . ,N is characterized by the sequence of bonds 
bn'-i, . . . , 6jv_i. Each bond represents a specific choice from the set of all possible 
bond vectors {b}. So, we can write the Rosenbluth weight of the monomers n', . . . ,N 
of chain b as 

W\x') = f[ w^{x') , w^{x') = e-'^-Ubi-i) + J2 e-'^"'^*^) . (51) 

Here, ?i* , (&„/-!) is the energy of monomer n' at its actual position in the chain. It in- 
cludes the interactions with the monomers of all other chains and with the monomers 
i = 1, . . . , n' — 1 of its own chain [Eq. (46)]. The monomers n' + lto N have to be 
omitted because one thinks of the chain b as being (re-) constructed step-by-step via 
the RR procedure with probability [Eq.(48)] 

w-[x ) 



^CBMC was introduced by Siepmann and Frenkel in the 1990s. It can be applied to lattice and off-lattice 
models. The initial off-lattice applications have demonstrated the power of the algorithm for the study of a large 
variety of problems in polymer physics. Therefore, CBMC has become an important and widely used simulation 
technique. A comprehensive and very pedagogical account of the method, including flowcharts of the algorithm 
and examples, is given in the textbook by Frenkel and Smit.'' 

^The RG algorithm was introduced in Refs. 104, 105. A detailed description may be found in Chap. 13.7 of the 
textbook by Frenkel and Smit.^' For practical appUcations, it is very helpful that the FORTRAN codes of the 
"Case Studies" may be downloaded from http : //molsim. chem.uva .nl/f renkel^mit. 
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Figure 16. Schematic of a CBMC move executed in a polymer solution. The initial configuration x' of the 
solution is shown in the left figure. From all chains of x' the chain b and one of its monomers, n', are chosen. 
Here, n' = N — 2. The chain portion n', . . . , A'^ is deleted (dashed grey lines connecting the open grey circles) 
and reconstructed step-by-step. Upon completion we obtain a new chain configuration, "chain a", and so a new 
configuration, x, of the solution. In the construction, the new bond 6i-i from monomer j — 1 to monomer 

1 {= N — 1, N) is selected from the set of bond vectors {b} according to the local Boltzmann factor P^^{x) 
[Eq. (55)]. For lattice models this set, which defines the number of trial directions at each step, is finite and closely 
related to lattice structure: E.g. for the hypercubic lattice, it coincides with the number of lattice directions z (or 

2 — 1 if backfolds shall be excluded a priori). For the BFM the trial directions may be taken as the ensemble of 
allowed bonds [Eq. (33)]. For a continuum model there is a priori an infinite number of possible directions, from 
which a suitable finite number of trial directions must be selected so that new chain configurations are efficiently 
generated (see Chap. 13.3 of Ref. 31 for details). If only excluded volume interactions are present, as assumed in 
the figure, P^^{x) = 1/2 for i = N — 1 and P^^{x) = 1/3 for the last monomer A'. Thus, the total weight 
W'^{x) of chain a is 6, whereas that of the chain portion n' = N — 1, N ot the old chain b is W^(x') = 4. 
Thus, W'^{x) /W^{x') = 1.5, and chain a will be accepted according to Eq. (58). 



Thus, the potential energy of the chain portion n' , . . . ,N is given by 

N 

U\x') = ^ u^(6,_i) , (53) 
and the probabihty to propose the chain portion may be expressed as 

N -f}U''{x') 



2. The monomers n' to N of chain h are deleted. This corresponds to a "shrinkage- 
growth" implementation of the algorithm, which is (presumably) more efficient than 
a "growth-shrinkage" procedure (see Sec. 5.2). 

3. To obtain a new configuration x of the system the chain h is fully (n! = 1) or only in 
part in' > 1) reconstructed: 

• If n' = 1, we start building a new chain by randomly placing the first monomer 
somewhere in system. Let this new chain be labeled "chain a". At the position, 
where the monomer is inserted, it interacts with the other chains of the system. 
It has an energy w5^=i(&o) [— U{fi), see Eq. (46)], giving rise to the weight 

wi ^ exp[-/3uJ(6o)]. 
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• If n' > 1, we rebuild the chain portion i = n', . . . ,N monomer-by-monomer 
via the RR method. This growing chain, comprising the initial portion of chain 
b (i ~ 1, . . . ,n' — 1) and the newly attached monomers (i = n, . . . , N) will also 
be called "chain a". In the following we suppress the prime (') to indicate that a 
new chain conformation is obtained, even if the first n—1 monomers are identical 
with those of chain b. To add the nth monomer we proceed as described in the 
next step. 

4. A monomer may be attached to the growing chain via a bond from the set of possible 
bond vectors {b}. Out of these trial directions we choose one according to Eq. (48). 
This impUes for the monomer i (= n, . . . , N) 

g-/3<(5i_i) g-/3«?(6i_i) 

^mI^') = = w^{x) ' ^^^^ 

where u°(6i_i) is the change in potential energy of the system due to the addition of 
the new bond from monomer i — 1 to monomer i [Eq. (46)]. 

5. The preceding step is repeated until the construction of the chain is completed. Thus, 
the new chain configuration occurs with probabiUty [Eq. (49)] 

i=n ^ ^ 

where W^^x) and U°'{x) are defined analogously to Eqs. (51,53). 

6. Now, we recognize that P^{x) may be interpreted as the probability of proposing a 
transition from the old configuration x' to the new configuration x. Correspondingly, 
P^{x') is the probability for the reverse step. That is, 

P^{x) - Pp„(a;' ^ x) and p!;{x') = Pp„(a; ^ x') . (57) 

Finally, we insert Eq. (57) into Eq. (29) to obtain the acceptance probability for the 
new configuration 

acc(a;' ^ a;) = min f 1, ^^"^"T ""'l e-^[^°(-)-^'(-')] 

Discussion. An important feature of the CBMC method is its non-local character: A suc- 
cessful CBMC construction implies a large-scale configurational change. Either a new 
chain is inserted somewhere in the system -this may be employed very efficiently to study 
phase equilibria of polymer solutions in the bulk^^* and in thin films^"^- or part^ of a chain 



^This part need not necessarily start at monomer n' + 1 and tenninate at the chain end, as assumed in our previous 
discussion. It can also comprise an inner portion of the chain, say the monomers n' + 1, . . . ,n' +m' — 1 < N. 
In this case, the reconstruction has to satisfy the additional condition that it starts at monomer n' and ends at the 
position of monomer n' + m' (see Chap. 13.4 of Ref. 31). This variant of CBMC may be used to relax e.g. ring 
polymers, which have no free ends so that the method described above could not be appUed. 
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is regrown. This leads to a rapid decorrelation of chain configurations and to efficient sam- 
pling, provided the system has a low or moderate density and the chains are not too long. 
When dealing with long chains and/or dense melts the following problems occur: 

• Chain construction in CBMC is based on the Rosenbluth method, yielding a distri- 
bution of configurations that differs from the Boltzmann distribution. This difference 
becomes more pronounced with increasing N, even for an isolated chain at the O- 
point (see Sees. 6.1,6.2). As CBMC does not control the weights while synthesizing 
the chain, contrary to the PERM, the acceptance rate for chain reconstructions faUs 
exponentially in N for large chain lengths. 

• Another problem results from the "shortsightedness" of the RR algorithm. By looking 
only one step ahead, the chain construction can run into traps. For isolated chains this 
trapping implies that there is still an exponential attrition in N for long chains, albeit 
with an attrition constant much smaller than A (see Eq. (44) and footnote on page 35). 

• In a dense system, in addition to trapping, a further problem occurs. If a (part of a) 
new chain is inserted, it is fairly likely to be constructed just in the space originally 
occupied by the (part of the) old chain which was removed. This can lead to strong 
correlations between the new and old chain configurations so that sampling becomes 
inefficient.^^ In this situation, it might be best to combine CBMC with the slithering- 
snake algorithm (Sec. 5.2). That is, to try to regrow just the terminal bond at the other 
chain end according to Eq. (58). In a dense melt, it thus appears as if CBMC cannot be 
expected to decorrelate chain configurations more efficiently than the slithering- snake 
algorithm. 

6.3.2 Recoil-Growth (RG) Algorithm 

The RG algorithm was suggested as an alternative to CBMC," exhibiting two major 
changes: 

1 . Instead of looking one step ahead the RG algorithm scans the environment via a re- 
tractable "feeler". The feeler consists of a self-avoiding chain portion having at most 
^recoil mouomcrs.'' The ability of the feeler to shrink and to grow helps to circum- 
vent dense regions. This allows for the search of suitable pathways to complete the 
construction of the chain. 



"The RG algorithm was introduced for SAW's on a cubic lattice in Ref. 104. Reference 105 extends this study 
to continuum models. A comprehensive discussion of these works may be found in the textbook by Frenkel and 
Smit (Chap. 13.7)?' 

''Of course, in general A^recoii > 2, implying that the feeler is longer than one bond. Otherwise, the "shortsight- 
edness" of CBMC is not removed. The idea to improve the RR method by looking several steps ahead is not 
new. It is embodied e.g. in the "scanning method" of Meirovitch. This method still uses a one-step growth, 
but chooses a new bond bi according to the probability that a SAW of Nsam monomers can be constructed in 
direction of 6j. As this impUes an enumeration of all possible SAW's of length iVscan starting at some monomer 
i, the scanning parameter ATscm is usually much smaller than N. Thus, trapping cannot be avoided completely. 
This would oiHy be the case if ATscan = N — i, i.e., if one scanned all possible ways to complete the chain up to 
monomer N in direction of 6j. For a comparative discussion of the RG algorithm and the scaiming method see 
Ref. 104. 
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2. Contrary to CBMC, the incremental weights iVj for each newly added monomer are 
not calculated "on the fly", but only after a new chain has been successfully con- 
structed. Thus, the computation of the weights is carried out only once. In CBMC, 
it can happen that a lot of time is spent to calculate the weights of a partially grown 
chain which must then be discarded because the construction has run into a trap before 
completion. 

Description of the RG Algorithm. These differences show that the RG algorithm comprises 
two independent steps: a construction step using the retractable feeler and an acceptance 
step including the weight calculation. We discuss these steps separately. In our discussion 
we assume that a whole chain is inserted in the system after a randomly chosen old chain 
has been removed ("shrinkage-growth procedure"). As before, the new chain is called 
"chain a", the old "chain 6", and the respective new and old configurations of the system 
are denoted by x and x'. 

Construction step (Fig. 17): 

1. We place the first monomer of chain a at the random trial position fi and determine 
its energy U{fi) = w?(&o)- To decide whether the position is accessible ("open") we 
accept it with probability (see e.g. Ref. 31 for further discussion of this point) 



where "SAW" means that there are only hard-core interactions.'^ If the position is 
"open", we continue with step 2. Otherwise, step 1 is repeated. 

2. Assume that the chain construction has arrived at monomer i. We randomly choose a 
new bond bi from monomer i to monomer compute (a;) for the trial position 
of i -h 1, and decide whether it is open or not according to Eq. (59). If it is closed, we 
keep choosing new bonds up to a maximum number of A;rg trial directions.'' As soon 
as the first open direction is found, we proceed to step 3. 

Otherwise, a recoil step is performed. That is, the chain moves back to monomer i — 1 
where it searches for an open direction, using the fcrest = fee — ^checked previously 
unchecked directions. The first open direction found is used to place (again) monomer 
i. If the chain fails to find an open direction from the fcrest remaining ones at monomer 
i — 1, it recoils to i — 2 and checks the previously unused directions for a possibility to 
grow from there. In difficult situations the chain keeps falling back until a maximum 
number of A^recoii recoil steps has been performed. If this number is exceeded, the 
construction of chain has to resume from step 1. 

''In the SAW-limit, the decision of whether the position is open is not probabilistic, but deterministic. One just 
checks on overlaps with other monomers: If no overlap occurs, the position is open. For continuous potentials, 
however, the decision becomes probabilistic. We compare P°^''{x) to a random number uniformly distributed 
between and 1. The position is open if C < P°''™(a;). Otherwise, it is closed. This means that we accept the 
position with probability P^'^{x). 

"^In practical implementations, the number of trial directions need not be the same for every monomer i. It may 
depend on i. For instance, one could choose for half of the monomers 2 trial directions and for the other half 3 
directions so that ^rg = 2.5 on average. Thus, feRo need not be integer. 




1 if fi is unoccupied , 
otherwise , 



(59) 
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Figure 17. Illustration of the recoil-growth procedure for a SAW on a squai'e lattice, using A;rg = 2 and Af,-ecoil = 
3. The first monomer is placed at an empty ("open") lattice site. At each site that the growing chain visits, it 
has fcRG trial directions for the next step (indicated by the dashed lines: there are 2 for each monomer). These 
directions are chosen at random from the z — 1 forward lattice directions. Starting at monomer 1 the chain first 
tries the path 1 2 3'. At that point, the first monomer (3 — A^rccoil + 1) becomes fixed (indicated by Q ^ » 'm the 
figure). Then, the chain attempts to grow to position 4' and finds it blocked. Thus, it recoils to monomer 3' and 
tries the remaining direction leading to monomer 4" . At that point, the monomer 4 — N^^con + 1 = 2 becomes 
fixed (O — > •). as the "feeler" 2 3'4" has attained its full length A^recoil = 3 (grey thick solid line). However, 
the next two trial directions, 4' 5' and 4' — ► 5", are found blocked. So, the feeler first recoils to 3', where 
it realizes that fcRo trial directions have been exhausted so that it must fall back its maximum length N^^co'd to 
monomer 2. From there, it finds the open path 2 3 4'", 5'" so that monomer 3 becomes permanently attached to 
the chain. In this way, the construction goes on up to monomer N — N^^^^a + 1. If this monomer was fixed, the 
construction stops, as a feeler to monomer A'^ exists. Figure adapted from Ref. 104. 



This shows that the RG algorithm is characterized by two parameters: k^c and A'^recoii- 
Tuning of these parameters is important to optimize the efficiency of the method. 

3. If the construction has successfully placed the nth monomer, we attach monomer 
n—A^,ecoii+l permanently to the chain, as recoiling can only fall backto n— A'recoii+l- 
If this occurs and still no open direction is found at monomer n — iVrecoii + 1, the 
growth process terminates and we must go back to step 1 . Thus, the A'recoii monomers, 
i — n — A'l-ecoii + 1, ■ ■ ■ ,n, may be regarded as "retractable feeler" of maximum length 
^recoil which probcs the territory ahead of monomer n 

-^recoil ■ 

4. Repeat steps 2 and 3 until the monomer N — A'lecoii + 1 was attached to the chain. 
This imphes that the feeler attained the chain end. A complete chain of length N has 
thus been constructed. 



Acceptance step: If chain a has been successfully constructed, we have to determine its 
weight W'^{x) and that of the old chain b, W''{x'), which chain a attempts to replace. 
This may be done as follows: 

1. In the construction step each monomer i disposes of (at most) /crg trial directions 
into which the growth of a feeler of maximum length A'lecoii may be attempted. As- 
sume that we have checked k ~ 1, . . . , /cchecked of these directions and found that the 
fcchecketjth direction is the first allowing us to complete the construction of the feeler 
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up to length A^recoii- We now test the remaining ktest = ^rg — ^checked directions to 
find out how many "feelers" of length A'recoii can be grown. (For the last monomers, 
N — A'recoii + 1 < i < A^, thc length of the feeler is shortened by one bond, as i 
advances step-by-step toward N.) For monomer i let mi(x) denote the number of 
directions where a feeler of full length can be grown. We have 1 < mj (x) < 1 -|- fcrest- 
As monomer i is only irrevocably added to the chain, if its position was initially open 
[Eq. (59)] and if at least one feeler of length iVrecoii may be grown from it, the proba- 
bility with which we propose to place monomer i is given by [see Eq. (55)] 

[x) = (^) • = -4^-^ ■ (60) 

With respect to CBMC P^^{x) differs by the absence of the Boltzmann factor 

exp[— /3w^(6i_i)], since only excluded- volume interactions are taken into account 
during the construction. So, the new chain a is proposed with probability 

« 1 1 
P°(a;) = TT — — = ——- iniNix) = 1) . (61) 

2. The weight calculation of the old chain b has two ingredients. First, we calculate 
P°'"^"(a;') for each monomer i according to Eq. (59).^ Second, we attempt to grow 
feelers in A:rg — 1 randomly chosen directions and we count the number of feelers that 
attain length A'recou- This number plus 1 yields mi{x') ("plus 1", since one feeler of 
length A'recou already exists along the backbone of chain 6). In analogy to chain a, we 
thus write 

^ pop™/ /\ N 

3. Finally, the total potential energy of chains a and h, [/"(a;) and lJ^{x'), must be 
computed before the new chain can be accepted with probability [see Eqs. (57,58)] 

acc(a=' - X) = min (l, ^-&\^^i^)-u\.-S\\^ . (gg) 

The RG algorithm has two adjustable parameters: ^rg, the (average) number of trial direc- 
tions of a monomer, and A^recoU, the length of the retractable feeler Intuitively, one expects 
that A^recoii should be large, whereas /jrg should be small. The main advantage of the RG al- 
gorithm is its ability to avoid traps by probing the environment with a feeler. A short feeler 
will strongly reduce this ability and thus the rate of successful chain construction. On the 
other hand, the value of /crg should not be too large because 1 < niiix) < k^c. (Remem- 
ber that m.i (x) denotes the number of feelers of length Ntecoih which may be grown from 
monomer i.) A large fcRc allows for many different values of mi{x), as we go along the 
chain to calculate the weight ^{x). The spread in mi{x) leads to a wide distribution of 
W°'{x), which will strongly deviate from the Boltzmarm distribution (see Sec. 6.1). 

'^For the old chain we only calculate P^^{x'). This is different from the constraction step, where we decide 
whether the position is open or closed according to probabiUty P°^^{x). That is, we first calculate P°''™(a;) 
and then compare it to a random number C, uniformly distributed between and 1. If C < P°''™(a;), the position 
is open. Otherwise, it is closed [see Eq. (59)]. 
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Figure 18. Left: Acceptance probability versus k^a for 8 < < 2048 obtained from simulations of isolated 
2D chains ("Kremer-Grest model") via the RG algorithm. For all cases, NrecoU = N/2. Note the sharp peak 
of acc(fcRG) for long chains. The value of fcRo at the peak yields the optimum choice fc^(W) for this chain 
length. For fe^^ we find: acc oc 1/N. Right: Plot of Rg /R^ versus N (main figure) and of R^ and ijg versus 
the number of bonds N — 1 (inset). The dashed horizontal Une in the main figure represents the large- AT limit of 
Rg /R^ [Eq. (64)]. The dotted line in the inset shows the theoretically predicted power law i?e oc JJg ~ AT^" 
with 1/ = 3/4. 



Applications of the RG Algorithm. These expectations are confirmed by applications of the 
RG algorithm to 3D lattice'""* and 3D off-lattice'"^ models of polymer solutions. For the 
chain lengths studied (A^ < 100) good choices are 2 < k^c < 3 and 3 < iVrecoii ^5 10' with 
the need to have larger values for both parameters if the density of the solution increases. 
Furthermore, it was observed that, for high densities and long chains, the RG algorithm 
may be an order of magnitude more efficient than CBMC. 

Our preliminary results'"^ on 2D polymer solutions, simulated with the Kremer-Grest 
model (see Sec. 4.2), appear to confirm the trends observed in 3D. Presumably due to the 
fact that the risk of trapping is more severe in 2D than in 3D,^^' ' '" we found that, even for 
an isolated chain, chain lengths TV > 100 are very difficult to simulate via CBMC. Here, 
the RG algorithm provides a powerful alternative. While the performance of the algorithm 
depends only weakly on the choice of A^recou. provided a long feeler is employed, e.g. 
-^recoil = N/2, fcRG must be optimized carefully. Figure 18 shows that the acceptance 
probability develops a pronounced peak for long chains. While for iV < 64 the acceptance 
probabiUty is fairly insensitive to the precise choice of fcRo, if /crg « 2, it has to be adjusted 
to 4 significant digits for N = 2048. In this case, the optimum value is = 1.838. This 
means that on average there are four monomers with two trial directions, followed by one 
monomer with only one trial direction. ^ 

To exempUfy that the RG algorithm produces correct statistical properties Fig. 18 also 
shows the radius of gyration R^, the end-to-end distance R^, and the ratio R^/R^. For 
isolated, two-dimensional chains one expects to find the critical exponent v = 3/4^^ and 



■'^We arrive at that conclusion in the following way. Suppose we allow the monomers to have either 1 or 2 
trial directions. Let p denote the probability that a monomer has two trial directions. Then, we pose feRo = 
2 • p H- 1 • (1 - p) so that p = 0.838 with fc"^! = 1.838 for N = 2048. Thus, out of 5 monomers, roughly four 
monomers have two trial directions, leaving one trial direction for the remaining monomer. 
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the ratio'" 

lim -I = 0.14026 ± 0.00011 . (64) 

The figure illustrates that both predictions are well borne out by the simulation data. 
Apparently, the 2D Kremer-Grest model quickly converges to the large-A^ limit, i.e, to 
i?e (X i?g ~ iV"^/^ and to Eq. (64). Deviations are only visible when investigating the 
ratio -RgZ-Rg for ^ ~ 512. For longer chains corrections-to-scaling are small. In this 
respect, the Kremer-Grest model agrees with the results obtained for the SAW on a square 
lattice."' 



7 Conclusions 

Behind the title of this chapter "Monte Carlo Simulation of Polymers: Coarse-Grained 
Models" a topic of a large breadth is hidden. So, a selection is necessary. We have em- 
ployed several criteria in this selection. 

First, we concentrate on simple generic models retaining only basic features of a poly- 
mer chain (chain connectivity, excluded volume, etc.; see Sec. 4). Coarse-grained models 
derived from specific polymers are only touched upon briefly (Appendix A), although this 
is an important current research topic.*'"''^ 

Second, we define a generic model as one consisting of monomers with the sim- 
plest imaginable structure. They are identified with sites on a lattice or with Lennard- 
Jones spheres in the continuum. The monomers in the chain are all the same and un- 
charged ("neutral homopolymers"). Their interaction is either purely repulsive or consists 
of a short-range repulsion supplemented by an attractive potential at larger distances (see 
Sec. 4). Thus, nor did we consider specific interactions, such as electrostatic interactions, 
H-bonds, interactions between different species of monomers (e.g. binary mixtures, block- 
copolymers), etc. -this will be done in other chapters of this book- neither did we dis- 
cuss current coarse-graining approaches which do not model a chain as a concatenation 
of monomers, but represent the whole polymer as a soft ellipsoidal"^'"^ or spherical"* 
particle. 

Within the scope of these generic models we presented various algorithms. What is the 
upshot of this discussion for applications? Here are some suggestions: 

• Dilute solutions: 

- For isolated chains, exempt of strong monomer-monomer interactions, the pivot 
algorithm is currently considered as being the most efficient method to study 
global properties related to the chain extension (iig, R^, the exponent v, etc.). 
To study properties related to the partition function (e.g. the exponent 7) other 
algorithms are better suited (see footnote on page 30). 

- To initialize the simulation via the pivot algorithm or for the analysis of local 
properties the slithering-snake algorithm represents an attractive alternative (see 
discussion on page 31). 

- If attractive monomer-monomer interactions are present, as it is the case close 
to the 6-point, it appears as if the pruned-enrichedRosenbluth method (PERM) 
is currently the most efficient algorithm (Sec. 6.2). 
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• From dilute solutions to dense melts: 



- As the density of the solution increases, the large-scale pivot moves become 
quickly inefficient. Here, either bilocal updating schemes, such as the sUthering- 
snake algorithm, or non-local chain reconstructions via configurational-bias 
Monte Carlo (CBMC) are better suited. CBMC has become a well established 
tool, particularly in its grand-canonical formulation for the study of phase dia- 
grams. It transpires that its range of applicability may be extended by the re- 
cently proposed recoil-growth (RG) algorithm (see Sec. 6.3). When the solution 
becomes more and more concentrated, the probability to renew the configura- 
tion of large chain portions becomes small. In this Umit, the CBMC and RG 
methods reduce to the slithering-snake algorithm. 

- In dense melts consisting of long chain (N > 10'^) also the shthering- snake 
algorithm faces serious problems to equiUbrate the system (Sec. 5.2). Here, it is 
better to use connectivity-altering moves instead of attempting to regrow (parts 
of) a chain. An example is the recently proposed double-pivot algorithm which 
appears to be very efficient in equilibrating dense melts (see Sec. 5.3). 

The previous points only refer to the study of conformational and structural properties 
of polymer systems as well as to an efficient equilibration of the system. If the focus 
of interest are dynamic properties, local moves should be employed because they mimic 
best the physical dynamics (if solvent-mediated hydrodynamic interactions are absent; see 
Ref. 14). 

The suggestions made above are a result of the discussion given in this chapter. Cer- 
tainly, our discussion suffers from omissions. Personally, we feel that the most serious 
one are generalized ensemble techniques,"^- such as simulated or parallel tempering 
(see Chap. 14.1 of Ref. 31 for an introduction). In the context of polymer physics, these 
methods have been applied e.g. to determine the chemical potential of polymer chains 
(simulated tempering)"^, to accelerate the equilibration of dense polymer melts (parallel 
tempering)"** or to the simulation of phase transitions"^. 

We apologize for this review-like end of our report, but hope that its content, together 
with the bibliography, wiU be helpful for those interested in MC simulations of polymer 
models. 

Appendix 

A Realistic Models and Coarse-Graining of Real Polymers 

The presentation of this chapter concentrated on generic coarse-grained models (Sec. 4). 
They are particularly useful for the study of the large-scale behavior of polymers. For these 
generic models several simulation algorithms were discussed. If one is now interested in 
properties of a specific molecule or material, the presented methods can still be applied. 
But their efficiency must be tested for the particular application.^^ The simulation of real 
materials must use models which reflect the underlying chemical structure of the polymer, 
even at the coarse-grained level. So, some mapping between a detailed, chemically realistic 
and a coarse-grained model must be established.'*^''^'^' ^-^^ The following paragraphs explain 
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the need and the main steps of such a coarse-graining procedure and briefly describe an 
application to modeling a melt of poly (vinyl alcohol). 

From Realistic to Coarse-Grained Models. An ideal simulation approach would consist 
in studying systems of long chains {N > 10^) with potentials being calculated from the 
simultaneous motion of all nuclei and electrons (via the Car-Parrinello method'^'). How- 
ever, already for an isolated chain such an approach is not feasible due to the large spread 
of subatomic and molecular time scales.''' ^ For A'' ~ 10^ the relaxation time of a polymer in 
dilute solution is about 1 |j.s, whereas the inclusion of the electrons requires a time step of 
~10~^ s in a MD simulation. This disparity of 1 1 orders of magnitude cannot be covered 
in present day simulations. A further restriction results from the system sizes that may 
be simulated. Typically, the total number of particles (electrons and nuclei) is limited to 
~103. 

Thus, simplifications are necessary.^* There are several levels to this.^^" In a first step, 
the "on the fly" calculated quantum-mechanical potentials can be replaced by empirical 
potentials for the bond lengths, the bond angles, the torsional angles and the non-bonded 
interactions between distant monomers along the backbone of the chain ("quantum level 

atomistic level"). After careful optimization of the parameters of these potentials quan- 
titative agreement between simulation and experimental results can be obtained.'*^' If 
the system is to stay in thermal equilibrium on both the local monomeric scale and the 
global scale of the chain, such comparisons between simulation and experiments are lim- 
ited to high temperatures up to now. Extensions to low temperatures where crystallization 
or sluggish glass-like relaxation occur still represent a great challenge for simulations on 
atomistic scale.'' To tackle these problems the models must be further simplified. In a 
second step of simplification, fast degrees of freedom (bond length and bond angle vibra- 
tions, etc.) may therefore be eliminated by a coarse-graining procedure ("atomistic level 

coarse-grained level"). This procedure leads either to generic models, as those which 
were discussed in the previous sections, or to a coarse-grained model for a specific polymer 
(see Sec. 4). Recent approaches to the latter case have been reviewed in Refs. 41, 42, 120. 

The first step, from the quantum to the atomistic level, may be interpreted as an ab-initio 
approach or a bottom-up construction of the model. When pushing the simplifications of 
the model beyond the atomistic to the coarse-grained level it becomes more and more 
apparent that the empirical potentials are just simulation parameters. So, instead of the 
bottom-up construction, one may also choose a more pragmatic top-down procedure. That 
is, given some (measured) properties, what are the potentials which can reproduce them? ' 
In practice, one often proceeds as follows: By the bottom-up method a first guess of coarse- 
grained potentials is constructed. In a second step, these potentials are optimized with 
respect to some experimentally or theoretically known properties. With the final potentials 

f Here, we assume that the whole chain is simulated at the same level of modeling. Of course, it is also possible to 
treat some part of the chain in atomistic detail (e.g. via Car-Parrinello), whereas others are modeled on a coarse- 
grained level. A recent example of such a multiscale approach is the selective adsorption of polycarbonate on a 
nickel surface. '^^ 

''An important contribution to this field are the "amorphous cell" simulations pioneered by Theodorou and 
Suter'^3' ™. Here, the idea is to fold an infinitely long chain in the simulation box such that the monomer density 
is close to the experimental one and the resulting chain structure is reasonable. With this approach mechanical 
properties of the amorphous, glassy melt have been studied successfully. 

*One such method for the solution of the inverse problem, known as "reverse Monte Carlo", is actually used to 
interpret neutron and x-ray scattering data.'^^''^* 
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Figure 19. Schematic illustration of the coarse-graining procedure developed in Ref. 127. Starting from an atom- 
istically detailed model of poly(vinyl alcohol) a bead-spring model is constructed. A monomer of the atomistic 
model is represented by a sphere on the coarse-grained level. The spheres interact by a harmonic bond poten- 
tial and a purely repulsive LJ-like potential [Eq. (65)]. The parameters of these potentials are optimized in the 
following way. In an atomistic simulation at high temperature, the structure of a melt of poly(vinyl alcohol) and 
the local conformation of a polymer are determined: the structure by the pair-distribution function g{r')''' and 
the conformation by the distribution of the angle between two vectors formed by connecting the C-atom of the 
CHOH-group of the (i — l)th monomer with that of the ith monomer (first vector) and that of the ith monomer 
with that of the (i + l)th monomer (second vector). From the angular distribution a potential of mean force" is 
calculated. In this case, it could directly be used as the bending potential for the bond angle in the coarse-grained 
model. Contrary to that, the parameters of the LJ-potentials are iteratively optimized until g(r) of the poly(vinyl 
alcohol) melt coincides with that of the coarse-grained model. 



one may then (try to) predict quantities which are not easily measurable. In the following 
we will sketch one such pragmatic approach. 

An Example of How to Coarse-Grain Real Polymers. The idea and the realization of the 
coarse-graining are illustrated by the example of poly(vinyl alcohol) in Fig. A.'^^ 

On the coarse-grained level, the chemical structure of the monomers is suppressed. 
They are represented by spheres, bound to each other by a harmonic potential. The orien- 
tation of successive bonds along the chain is determined by a bond angle potential, whereas 
other non-bonded monomers interact by the repulsive part of a LJ-like potential 



where r^i^ is the minimum of the potential and eo is chosen such that the minimum is 
shifted to zero to have a continuous force. When fitting the potential to real polymers it 




min J 



(65) 



50 



turned our^ that for the repulsive term, an exponent smaller than 12 may be better adapted 
because the coarse-grained units have to be somewhat softer than atoms.^ By comparing 
high-temperature simulations of an atomistically detailed and a coarse-grained model, the 
potential parameters of the latter are adjusted in such a way that it faithfully reproduces 
the backbone rigidity of poly(vinyl alcohol) and the local packing of its monomers in the 
liquid state. 

The construction of the model was done at one state point. It is a priori not guaranteed 
how far from this reference conditions the model still matches the behavior of the original 
one. It turns out that in the present case, the crystallization at lower temperatures is well 
reproduced. 
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